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Preface 



These notes were originally prepared during Fall quarter 2003 for UCSD Math 174, Numerical 
Methods. In writing these notes, it was not my intention to add to the glut of Numerical Analysis 
texts; they were designed to complement the course text, Numerical Mathematics and Computing, 
Fourth edition, by Cheney and Kincaid [7]. As such, these notes follow the conventions of that 
text fairly closely If you are at all serious about pursuing study of Numerical Analysis, you should 
consider acquiring that text, or any one of a number of other fine texts by e.g., Epperson, Hamming, 
etc. [3, 4, 5]. 



Figure 1: The chapter dependency of this text, though some dependencies are weak. 

Special thanks go to the students of Math 174, 2003-2004, who suffered through early versions 
of these notes, which were riddled with (more) errors. 

Revision History 

0.0 Transcription of course notes for Math 174, Fall 2003. 
0.1 As used in Math 174, Fall 2004. 
0.11 Added material on functional analysis and Orthogonal Least Squares. 
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• More homework questions and example problems. 

• Chapter on optimization. 

• Chapters on basic finite difference and finite element methods? 

• Section on root finding for functions of more than one variable. 
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Chapter 1 

Introduction 



1.1 Taylor's Theorem 

Recall from calculus the Taylor's series for a function, f(x), expanded about some number, c, is 
written as 

f(x) ~ a + ai (x - c) + a 2 (x - c) 2 + 

Here the symbol ~ is used to denote a "formal series," meaning that convergence is not guaranteed 
in general. The constants a,{ are related to the function / and its derivatives evaluated at c. When 
c = 0, this is a MacLaurin series. 

For example we have the following Taylor's series (with c = 0): 



2 S 

sin (x) = x - — + — -.. . 

x 2 x 4 
cos(x) = l-- + ir -... 



(1.1) 

(1.2) 
(1.3) 



Theorem 1.1 (Taylor's Theorem). If /(x) has derivatives of order 0, 1, 2, . . . , n+1 on the closed 
interval [a, 6] , then for any x and c in this interval 



f(x) = f(k) (c) {x ~ c)k + /(n+1) ^ {x ~ c)n+1 



k=0 



k\ 



(n + 1)! 



where ^ is some number between x and c, and f k (x) is the /c th derivative of / at x. 



We will use this theorem again and again in this class. The main usage is to approximate 
a function by the first few terms of its Taylor's series expansion; the theorem then tells us the 
approximation is "as good" as the final term, also known as the error term. That is, we can make 
the following manipulation: 
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k=0 



/W (c) (x-c)< 



fc=0 



fc! 



/W ( 



C ) {X-C) k (0 (X-C) 

fc! 



E 

fc=0 

/(n+D (g) (g _ c) n+l 

(n + 1)! 

1/^) (oiix- c r +i 

(n + 1)! 



n+1 



(n+1)! 



On the left hand side is the difference between f(x) and its approximation by Taylor's series. 
We will then use our knowledge about /( n+1 ) (£) on the interval [a, b] to find some constant M such 
that 



/(*) - E 



c x — c 



fc=0 



k\ 



|/(n+D (Q| |x-C 

(n + 1)! 



n+1 



< M |x - c 



n+1 



Example Problem 1.2. Find an approximation for f(x) = sinx, expanded about c = 0, using 
n = 3. 

Solution: Solving for is fairly easy for this function. We find that 



f(x) = sinx = sin(O) + 



cos(O) x — sin(0) x 2 — cos(O) x 3 sin(£) x 



1! 



+ 



2! 



+ 



3! 



+ 



4! 



so 



x 3 sin(f ) x 4 



sinx 



sin(£) x 4 



24 



x 4 
< — , 
- 24' 



because |sin(£)| < 1. 



Example Problem 1.3. Apply Taylor's Theorem for the case n = 1. 

Solution: Taylor's Theorem for n = 1 states: Given a function, /(x) with a continuous derivative 
on [a,b], then 

m = f(c)+f(o(x-c) 

for some £ between x, c when x, c are in [a, b]. 

This is the Mean Value Theorem. As a one-liner, the MVT says that at some time during a trip, 
your velocity is the same as your average velocity for the trip. H 



Example Problem 1.4. Apply Taylor's Theorem to expand f(x) 
Solution: Simple calculus gives us 



21x + 17 around c = 1. 



\x 
'(x 
'(x 
)(x 
)(x 



= x 3 - 21x 2 + 17, 

= 3x 2 - 42x, 

= 6x - 42, 

= 6, 

= 0. 
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with the last holding for k > 3. Evaluating these at c = 1 gives 
f(x) = -3 + -39(x - 1) + y - '- + 



6 

Note there is no error term, since the higher order derivatives are identically zero. By carrying out 
simple algebra, you will find that the above expansion is, in fact, the function f(x). H 

There is an alternative form of Taylor's Theorem, in this case substituting x + h for x, and x 
for c in the more general version. This gives 

Theorem 1.5 (Taylor's Theorem, Alternative Form). If f(x) has derivatives of order 0, 1, ... , n+ 
1 on the closed interval [a,b], then for any x in this interval and any h such that x + h is in this 
interval, 

/(l + ft) ^Z!^ + /''7«)(f +1 , 

^ fc! (n+1)! 

where £ is some number between x and x + h. 

We generally apply this form of the theorem with h — > 0. This leads to a discussion on the 
matter of Orders of Convergence. The following definition will suffice for this class 

Definition 1.6. We say that a function f(h) is in the class O (h k ) (pronounced "big-Oh of h k ") 
if there is some constant C such that 

\f(h)\<C\h\ k 

for all h "sufficiently small," i.e., smaller than some h* in absolute value. 

For a function / G O (h k ) we sometimes write / = O (/i fc ) . We sometimes also write O (h h ), 
meaning some function which is a member of this class. 

Roughly speaking, through use of the "Big-O" function we can write an expression without 
"sweating the small stuff." This can give us an intuitive understanding of how an approximation 
works, without losing too many of the details. 

Example 1.7. Consider the Taylor expansion of Inx: 

w u . , (i/x) h (-i/* 2 ) /i 2 (2/e) 

\n(x + h) = lnx+ v/ i 7 + ^ '-—I + 

Letting x = 1, we have 

ln(l + h) = h-^ + ^h\ 

Using the fact that £ is between 1 and 1 + h, as long as h is relatively small (say smaller than ^), 
the term can be bounded by a constant, and thus 

h 2 

In (1 + ft) = h- — + 0(h 3 ). 



Thus we say that h — 4^ is a O (/i 3 ) approximation to ln(l + h). For example 

01 2 

ln(l + 0.01) 0.009950331 w 0.00995 = 0.01 . 
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1.2 Loss of Significance 

Generally speaking, a computer stores a number x as a mantissa and exponent, that is x = ±r x 10 fc , 
where r is a rational number of a given number of digits in [0.1, 1), and k is an integer in a certain 
range. 

The number of significant digits in r is usually determined by the user's input. Operations 
on numbers stored in this way follow a "lowest common denominator" type of rule, i.e., precision 
cannot be gained but can be lost. Thus for example if you add the two quantities 0.171717 and 
0.51, then the result should only have two significant digits; the precision of the first measurement 
is lost in the uncertainty of the second. 

This is as it should be. However, a loss of significance can be incurred if two nearly equal 
quantities are subtracted from one another. Thus if I were to direct my computer to subtract 
0.177241 from 0.177589, the result would be .348 x 10~ 3 , and three significant digits have been lost. 
This loss is called subtractive cancellation, and can often be avoided by rewriting the expression. 
This will be made clearer by the examples below. 

Errors can also occur when quantities of radically different magnitudes are summed. For exam- 
ple 0.1234 + 5.6789 x 10 -20 might be rounded to 0.1234 by a system that keeps only 16 significant 
digits. This may lead to unexpected results. 

The usual strategies for rewriting subtractive expressions are completing the square, factoring, 
or using the Taylor expansions, as the following examples illustrate. 

Example Problem 1.8. Consider the stability of \Jx + 1 — 1 when x is near 0. Rewrite the 
expression to rid it of subtractive cancellation. 

Solution: Suppose that x = 1.2345678 x 10~ 5 . Then ^/xTT ~ 1.000006173. If your computer 
(or calculator) can only keep 8 significant digits, this will be rounded to 1.0000062. When 1 is 
subtracted, the result is 6.2 x 10 -6 . Thus 6 significant digits have been lost from the original. 
To fix this, we rationalize the expression 



r~ T a i — XT .v^TT+l x + 1-1 
\Jx + 1 — 1 = yx + l — l- 



'Vs+T+l Vx + 1 + 1 Vx + 1 + 1 ' 

This expression has no subtractions, and so is not subject to subtractive cancelling. When x = 
1.2345678 x 10 -5 , this expression evaluates approximately as 

1.2345678 x 10~ 5 fi 

» 6.17281995 x 10~ 6 

2.0000062 

on a machine with 8 digits, there is no loss of precision. H 

Note that nearly all modern computers and calculators store intermediate results of calculations 
in higher precision formats. This minimizes, but does not eliminate, problems like those of the 
previous example problem. 

Example Problem 1.9. Write stable code to find the roots of the equation x 2 + bx + c = 0. 
Solution: The usual quadratic formula is 

-b ± yjb 2 - 4c 

X ± = o 



Supposing that & > c> 0, the expression in the square root might be rounded to b 2 , giving two 
roots x + = 0, X- = —b. The latter root is nearly correct, while the former has no correct digits. 



1.3. VECTOR SPACES, INNER PRODUCTS, NORMS 



5 



To correct this problem, multiply the numerator and denominator of x + by —b — \Jb 2 — 4c to get 

2c 



x + 



-b - Vb 2 - 4c 



Now if b » c > 0, this expression gives root x+ = —c/b, which is nearly correct. This leads to the 
pair: 

-b - Vb' 2 - 4c 2c 

o > X + 



X_ = 



-6 - V6 2 - 4c 



Note that the two roots are nearly reciprocals, and if x- is computed, x + can easily be computed 
with little additional work. H 

Example Problem 1.10. Rewrite e x — cosx to be stable when x is near 0. 
Solution: Look at the Taylor's Series expansion for these functions: 



cosx 



tZ/^ tZJ^ 

1 + X+ ¥ + 3! + ¥ + 5f + 



= x + x 2 + ||- + C (x 5 ) 

This expression has no subtractions, and so is not subject to subtractive cancelling. Note that we 
propose calculating x + x 2 + x 3 /6 as an approximation of e x — cosx, which we cannot calculate 
exactly anyway. Since we assume x is nearly zero, the approximate should be good. If x is very 
close to zero, we may only have to take the first one or two terms. If x is not so close to zero, we 
may need to take all three terms, or even more terms of the expansion; if x is far from zero we 
should use some other technique. H 



2 4 fi 

X X X 

~ 21 + IT ~ 6f + 



1.3 Vector Spaces, Inner Products, Norms 

We explore some of the basics of functional analysis which may be useful in this text. 
1.3.1 Vector Space 

A vector space is a collection of objects together with a binary operator which is defined over an 
algebraic field. 1 The binary operator allows transparent algebraic manipulation of vectors. 

Definition 1.11. A collection of vectors, V, with a binary addition operator, +, defined over V, 
and a scalar multiply over the real field R, forms a vector space if 

1. For each u, v G V, the sum u + v is a vector in V. (i.e., the space is "closed under addition.") 

2. Addition is commutative: u + v = v + u for each u, v G V. 

3. For each u G V, and each scalar a G R the scalar product au is a vector in V. (i.e., the space 
is "closed under scalar multiplication.") 

4. There is a zero vector G V such that for any u G V, Ou = 0, where is the zero of R. 

5. For any u G V, lu = u, where 1 is the multiplicative identity of R. 

6. For any u, v G V, and scalars a, (3 G R, both (a +k [3)u = au + [3u and a (u + v) = au + av 
hold, where +m is addition in R. (i.e., scalar multiplication distributes in both ways.) 



1 For the purposes of this text, this algebraic field will always be the real field, E, though in the real world, the 
complex field, C, has some currency. 
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Example 1.12. The most familiar example of a vector space is R ra , which is the collection of 
n-tuples of real numbers. That is u G R n is of the form [ui, u 2 , ■ ■ ■ , u n ] T , where iij £ 1 for 
i = 1, 2, . . . , n. Addition and scalar multiplication over R are defined pairwise: 

[ui,u 2 ,.. • ,n n ] T + [vi,v 2 , . . . ,v n ] T = [ui + v 1 ,u 2 + v 2 ,.. -,u n + v n ] T , and 
a [ui,u 2 , . . . ,u n ] T = [aui,av,2, ■ ■ ■ ,au n ] T 

Note that some authors distinguish between points in n-dimensional space and vectors in re- 
dimensional space. We will use R n to refer to both of them, as in this text there is no need to 
distinguish them symbolically. 

Example 1.13. Let X C R fc be an closed, bounded set, and let H be the collection of all functions 
from X to R. Then H forms a vector space under the "pointwise" defined addition and scalar 
multiplication over R. That is, for u, v G H, u + v is the function in H defined by [u + v] (x) = 
u(x) + v(x) for all x G X. And for u G H, a G R, era is the function in [era] (x) = a (u(x)). 

Example 1.14. Let X C R fc be a closed, bounded set, and let Hq be the collection of all functions 
from X to R that take the value zero on dX. Then H forms a vector space under the "pointwise" 
defined addition and scalar multiplication over R. The only difference between proving Hq is a 
vector space and the proof required for the previous example is in showing that Hq is indeed closed 
under addition and scalar multiplication. This is simple because if x G dX, then [u + v] (x) = 
u(x) + v(x) = + 0, and thus u + v has the property that it takes value on dX. Similarly for au. 
This would not have worked if the functions of Hq were required to take some other value on dX, 
like, say, 2 instead of 0. 

Example 1.15. Let V n be the collection of all 'formal' polynomials of degree less than or equal to 
n with coefficients from R. Then V n forms a vector space over R. 

Example 1.16. The collection of all real-valued mxn matrices forms a vector space over the reals 
with the usual scalar multiplication and matrix addition. This space is denoted as R mxn . Another 
way of viewing this space is it is the space of linear functions which carry vectors of R n to vectors 
of R m . 

1.3.2 Inner Products 

An inner product is a way of "multiplying" two vectors from a vector space together to get a scalar 
from the same field the space is defined over (e.g., a real or a complex number). The inner product 
should have the following properties: 

Definition 1.17. For a vector space, V, defined over R, a binary function, (,), which takes two 
vectors of V to R is an inner product if 

1. It is symmetric: (v,u) = (u,v). 

2. It is linear in both its arguments: 

(au + [3v, w) = a (u, w) + [3 (v, w) and 
(u, av + f3w) = a (u, v) + (3 (u, w) . 

A binary function for which this holds is sometimes called a bilinear form. 

3. It is positive: (v, v) > 0, with equality holding if and only if v is the zero vector of V. 
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Example 1.18. The most familiar example of an inner product is the L 2 (pronounced "L two") 
inner product on the vector space R n . If u = [u±, 112, ■ ■ ■ , u n ] T , and v = [v 1, V2, ■ ■ ■ , v n ] T , then 
letting 

(u,v) 2 = y^ujVj 

i 

gives an inner product. This inner product is the usual vector calculus dot product and is sometimes 
written as u ■ v or u v. 

Example 1.19. Let H be the vector space of functions from X to R from Example 1.13. Then 
for u,v G H, letting 

{u,v) H = I u(x)v(x)dx, 
Jx 

gives an inner product. This inner product is like the "limit case" of the L 2 inner product on R n 
as n goes to infinity. 

1.3.3 Norms 

A norm is a way of measuring the "length" of a vector: 

Definition 1.20. A function ||-|| from a vector space, V, to R + is called a norm if 

1. It obeys the triangle inequality: \\x + y\\ < \\x\\ + \\y\\ . 

2. It scales positively: ||aa;|| = |a| , for scalar a. 

3. It is positive: > 0, with equality only holding when x is the zero vector. 

The easiest way of constructing a norm is on top of an inner product. If (, ) is an inner product 
on vector space V, then letting 

||ti|| = yj(u,u) 

gives a norm on V. This is how we construct our most common norms: 
Example 1.21. For vector x G R™, its L 2 norm is defined 



n \ 2 



This is constructed on top of the L 2 inner product. 



Example 1.22. The L p norm on R n generalizes the L 2 norm, and is defined, for p > 0, as 

/ n \ 1 /P 



X 



1= EN' 



\i=l 

Example 1.23. The L°° norm on R n is defined as 



\ x \\™ = max . 

00 i 



The L°° norm is somehow the "limit" of the L p norm as p —>■ 00. 
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1.4 Eigenvalues 

It is assumed the reader has some familiarity with linear algebra. We review the topic of eigenvalues. 

Definition 1.24. A nonzero vector x is an eigenvector of a given matrix A, with corresponding 
eigenvalue A if 

Ax = Xx 

Subtracting the right hand side from the left and gathering terms gives 

(A — X\)x = 0. 

Since x is assumed to be nonzero, the matrix A — Al must be singular. A matrix is singular if and 
only if its determinant is zero. These steps are reversible, thus we claim A is an eigenvalue if and 
only if 

det (A - Al) = 0. 

The left hand side can be expanded to a polynomial in A, of degree n where A is an n x n matrix. This 
gives the so-called characteristic equation. Sometimes eigenvectors, -values are called characteristic 
vectors,-values. 



Example Problem 1.25. Find the eigenvalues of 

" 1 1 

4 -2 



Solution: The eigenvalues are roots of 



= det 



1 - A 
4 



1 



-2- A 

This equation has roots Ai = — 3, A2 = 2. 



(l-A)(-2-A)-4 = A 2 + A-6. 



Example Problem 1.26. Find the eigenvalues of A 2 . 

Solution: Let A be an eigenvalue of A, with corresponding eigenvector x. Then 

A 2 x = A (Ax) = A (Ax) = AAa; = A 2 cc. 



The eigenvalues of a matrix tell us, roughly, how the linear transform scales a given matrix; the 
eigenvectors tell us which directions are "purely scaled." This will make more sense when we talk 
about norms of vectors and matrices. 



1.4.1 Matrix Norms 

Given a norm ||-|| on the vector space, M. n , we can define the matrix norm "subordinate" to it, as 
follows: 

Definition 1.27. Given a norm ||-|| on M. n , we define the subordinate matrix norm on IR nxn by 

IIAxll 
A = max —r — TT-. 

x^O \\X\\ 
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We will use the subordinate two-norm for matrices. From the definition of the subordinate 
norm as a max, we conclude that if a; is a nonzero vector then 



|Ax| 



- < ||A|| 2 thus, 



| J\tjC 1 1 2 1 1 1 1 2 1 1 *^ 1 1 2 * 



Example 1.28. Strange but true: If A is the set of eigenvalues of A, then 



Example Problem 1.29. Prove that 



I AIL = max I A| . 

uz AeA 



|AB|| 2 < ||A|| 2 ||B|| 2 . 



Solution: 

• I I I A B X cy /\ IL II BX|L Mam i i i i 

AB L =df max — — — < max - — ^- Tl — — = A L B L. 

x^O X 2 x^O X L 
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Exercises 



(1.1) 



1.2 
1.3 
1.4 
1.5 
1.6 
1.7 
1.8 
1.9 



1.10) 
1.11) 
1.12) 
1.13) 

1.14) 

1.15) 
1.16) 

1.17) 

1.18) 

1.19) 

1.20) 

1.21) 
1.22) 
1.23) 
1.24) 



Suppose / G O (h k ) . Show that / G O (h m ) for any m with < m < k. (Hint: Take 
h* < 1.) Note this may appear counterintuitive, unless you remember that O (h k ) is a better 
approximation than O (h m ) when m < k. 

Suppose f £ O (h k ) , and g G O (ft™) . Show that fg £ O (h k+m ) . 

Suppose / eO (/i fc ) , and g £ O (h m ) , with m < fc. Show that f + g £ O (h m ) . 

Prove that /(/i) = -3/i 5 is in O (h 5 ). 

Prove that f(h) = h 2 + 5/i 17 is in (/i 2 ). 

Prove that f(h) = h 3 is not in (/i 4 ) (Hint: Proof by contradiction.) 

Prove that sin(/i) is in (h). 

Find a (h 3 ) approximation to sin/i. 

Find a (/i 4 ) approximation to ln(l + /i). Compare the approximate value to the actual when 
h = 0.1. How does this approximation compare to the O (/j 3 ) approximate from Example 1.7 

for h = o.n 

Suppose that / G O (h k ). Can you show that /' G O (h k ~ l )l 

\/l to get rid of subtractive cancellation when x « 0. 
yfx to get rid of subtractive cancellation when x is very large. 

cosx of subtractive cancellation for x 

5 ^ approximate. 

cos 2 x of subtractive cancellation for x 



Rewrite \fx + 1 
Rewrite \Jx + 1 

Use a Taylor's expansion to rid the expression 1 
small. Use a O (x 5 ) approximate. 

Use a Taylor's expansion to rid the expression 1 - 
small. Use a (x 6 ) approximate. 
Calculate cos(7r/2 + 0.001) to within 8 decimal places by using the Taylor's expansion. 
Prove that if x is an eigenvector of A then ax is also an eigenvector of A, for the same 
eigenvalue. Here a is a nonzero real number. 

Prove, by induction, that if A is an eigenvalue of A then A fc is an eigenvalue of A fc for integer 
k > 1. The base case was done in Example Problem 1.26. 

Let B = Yli=o ai ^ l > wnere A = I. Prove that if A is an eigenvalue of A, then X^=o ai/ ^ ^ s 
an eigenvalue of B. Thus for polynomial p(x), p(X) is an eigenvalue of p(A). 

Suppose A is an invertible matrix with eigenvalue A. Prove that A -1 is an eigenvalue for 
A" 1 . 

Suppose that the eigenvalues of A are 1, 10, 100. Give the eigenvalues of B = 3A 3 — 4A 2 + I. 
Show that B is singular. 
Show that if || x \\ 2 = r, then x is on a sphere centered at the origin of radius r, in M. n . 
If ||x|| 2 = 0, what does this say about vector xl 



Letting x = [3 4 12] 1 
What is the norm of 



what is IIjcIL? 



A = 



1 







1/2 










1/3 









1/n 



(1.25) Show that ||A|| 2 = implies that A is the matrix of all zeros. 

(1.26) Showthat || A" 1 1| 2 equals (1/| A 

mini) > where \ m in is the smallest, in absolute value, eigenvalue 

of A. 

(1.27) Suppose there is some \i > such that, for a given A, 



|At?|| 2 > fj,\\v\ 



2- 
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for all vectors v. 

(a) Show that fx < ||A|| 2 . (Should be very simple.) 

(b) Show that A is nonsingular. (Recall: A is singular if there is some x / such that 
Ax = 0.) 

(c) Show that ||A _1 || 2 < (l//x) . 

(1.28) If A is singular, is it necessarily the case that ||A|| 2 = 0? 

(1.29) If ||A|| 2 > /i > does it follow that A is nonsingular? 

(1.30) Towards proving the equality in Example Problem 1.28, prove that if A is the set of eigen- 
values of A, then 

IIAII > max |A| , 

AeA 

where ||-|| is any subordinate matrix norm. The inequality in the other direction holds when 
the norm is ||-|| 2 , but is difficult to prove. 
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Chapter 2 

A "Crash" Course in octave/Matlab 



2.1 Getting Started 

Matlab is a software package that allows you to program the mathematics of an algorithm without 
getting too bogged down in the details of data structures, pointers, and reinventing the wheel. 
It also includes graphing capabilities, and there are numerous packages available for all kinds of 
functions, enabling relatively high-level programming. Unfortunately, it also costs quite a bit of 
money, which is why I recommend the free Matlab clone, octave, available under the GPL 1 , freely 
downloadable from http : / /www . octave . org. 

In a lame attempt to retain market share, Mathworks continues to tinker with Matlab to make 
it noncompatible with octave; this has the side effect of obsoletizing old Matlab code. I will try 
to focus on the intersection of the two systems, except where explicitly noted otherwise. What 
follows, then, is an introduction to octave; Matlab users will have to make some changes. 

You can find a number of octave/Matlab tutorials for free on the web; many of them are 
certainly better than this one. A brief web search reveals the following excellent tutorials: 

• http : //www . math .mtu . edu/~msgocken/intro/ intro . html 

• http : //www. cyclismo . org/tutorial/matlab/vector .html 

• http : //web . ew . usna . edu/~mecheng/DESIGN/CAD/MATLAB/usna . html 

Matlab has some demo programs covering a number of topics- from the most basic functionality 
to the more arcane toolboxes. In Matlab, simply type demo. 

What follows is a lame demo for octave. Start up octave. You should get something like: 

GNU Octave, version 2.1.44 (i686-pc-linux-gnu) . 

Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003 John W. Eaton. 
This is free software; see the source code for copying conditions. 
There is ABSOLUTELY NO WARRANTY; not even for MERCH ANT I B I L I TY or 
FITNESS FOR A PARTICULAR PURPOSE. For details, type 'warranty'. 

Please contribute if you find this software useful. 

For more information, visit http://www.octave.org/help-wanted.html 

Report bugs to <bug-octave@bevo . che . wise . edu> . 

octave : 1> 

1 Gmi Public License. See http://www.gnu.org/copyleft/gpl.html. 
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You now have a command line. The basic octavian data structure is a matrix; a scalar is a 1 x 1 
matrix, a vector is an n x 1 matrix. Some simple matrix constructions are as follows: 

octave :2> a = [12 3] 
a = 

12 3 

octave :3> b = [5; 4; 3] 
b = 

5 

4 

3 

octave :4> c = a' 
c = 

1 

2 

3 

octave :5> d = 5*c - 2 * b 
d = 
-5 

2 

9 

You should notice that octave "echoes" the lvalues it creates. This is either a feature or an 
annoyance. It can be prevented by appending a semicolon at the end of a command. Thus the 
previous becomes 

octave :5> d = 5*c - 2 * b; 
octave : 6> 

For illustration purposes, I am leaving the semicolon off. To access an entry or entries of a 
matrix, use parentheses. In the case where the variable is a vector, you only need give a single 
index, as shown below; when the variable is a matrix, you need give both indices. You can also 
give a range of indices, as in what follows. 

WARNING: vectors and matrices in octave/Matlab are indexed starting from 1, and not 
from 0, as is more common in modern programming languages. You are warned! Moreover, the 
last index of a vector is denoted by the special symbol end. 

octave :6> a(l) = 77 
a = 

77 2 3 

octave :7> a(end) = -400 
a = 

77 2 -400 
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octave :8> a(2:3) = [22 333] 



a = 



77 22 333 

octave :9> M = diag(a) 
M = 



77 











22 











333 



octave :10> M(2,l) = 14 
M = 

77 
14 22 
333 

octave:ll> M(l:2,l:2) = [1 2;3 4] 
M = 

12 
3 4 
333 

The command diag(v) returns a matrix with v as the diagonal, if v is a vector. diag(M) returns 
the diagonal of matrix M as a vector. 

The form c : d gives returns a row vector of the integers between a and d, as we will examine 
later. First we look at matrix manipulation: 

octave: 12> j = M * b 
j = 



13 
31 
999 



octave : 13> 



N = rand (3, 3) 



N = 



0.706307 
0.911833 



0.166880 



0.027866 0.087402 
0.624716 0.067067 
0.769423 0.938714 



octave : 14> 
L = 



L = 



M + N 



1.166880 
3.706307 
0.911833 



2.027866 0.087402 
4.624716 0.067067 
0.769423 333.938714 



octave :15> P = L * M 
P = 
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7.2505e+00 
1.7580e+01 
3.2201e+00 



1.0445e+01 
2.5911e+01 
4.9014e+00 



2.9105e+01 
2.2333e+01 
1.1120e+05 



octave :16> P = 
P = 

1.1669e+00 
1.1119e+01 
0.0000e+00 



L .* M 

4.0557e+00 
1.8499e+01 
0.0000e+00 



0.0000e+00 
0.0000e+00 
1.1120e+05 



octave : 17> x = M \ b 
x = 

-6.0000000 

5.5000000 

0.0090090 



octave :18> err = M * x - b 
err = 







Note the difference between L * M and L . * M; the former is matrix multiplication, the latter 
is element by element multiplication, i.e., 

(L-*M). j = H j K id . 

The command rand(m,n) gives an m x n matrix with each element "uniformly distributed" on 
[0, 1]. For a zero mean normal distribution with unit variance, use randn(m,n) . 
In line 17 we asked octave to solve the linear system 

Mx = b, 

by setting 

x = M\b = M _1 b. 

Note that you can construct matrices directly as you did vectors: 

octave :19> B = [1 3 4 5; 2 -2 2 -2] 
B = 

13 4 5 
2-2 2-2 

You can also create row vectors as a sequence, either using the form c : d or the form c : e : d, which 
give, respectively, c, c + 1, . . . , d, and c, c + e, . . . , d, (or something like it if e does not divide d — c) 
as follows: 



octave :20> z = 1:5 
z = 



1 2 3 4 5 
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octave:21> z = 5: (-1) : 1 
z = 

5 4 3 2 1 

octave :22> z = 5:(-2):l 
z = 

5 3 1 

octave :23> z = 2:3:11 
z = 

2 5 8 11 

octave :24> z = 2:3:10 
z = 

2 5 8 

Matrices and vectors can be constructed "blockwise." Blocks in the same row are separated by a 
comma, those in the same column by a semicolon. Thus 

octave :2> y=[2 7 9] 

y = 

2 7 9 

octave :3> m = [z;y] 
m = 

2 5 8 

2 7 9 

octave :4> k = [(3:4) ' , m] 
k = 

3 2 5 8 

4 2 7 9 



2.2 Useful Commands 

Here's a none too complete listing of useful commands in octave: 

• help is the most useful command. 

• floor (X) returns the largest integer not greater than X. If X is a vector or matrix, it computes 
the floor element-wise. This behavior is common in octave: many functions which we normally 
think of as applicable to scalars can be applied to matrices, with the result computed element- 
wise. 

• ceil(X) returns the smallest integer not less than X, computed element-wise. 

• sin(X), cos(X), tan(X), atan(X), sqrt(X), returns the sine, cosine, tangent, arctangent, 
square root of X, computed elementwise. 

• exp(X) returns e x , elementwise. 

• abs(X) returns |X| , elementwise. 
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norm(X) returns the norm of X; if X is a vector, this is the L 2 norm: 

/ \ 1/2 



if X is a matrix, it is the matrix norm subordinate to the L 2 norm. 

You can compute other norms with norm(X,p) where p is a number, to get the L p norm, or 
with p one of Inf , -Inf , etc. 

zeros (m,n) returns an m x n matrix of all zeros, 
eye (m) returns the m x m identity matrix. 

[m,n] = size (A) returns the number of rows, columns of A. Similarly the functions rows (A) 

and columns (A) return the number of rows and columns, respectively. 

length (v) returns the length of vector v, or the larger dimension if v is a matrix. 

f ind(M) returns the indices of the nonzero elements of N. This may not seem helpful at first, 

but it can be very useful for selecting subsets of data because the indices can be used for 

selection. Thus, for example, in this code 

octave :1> v = round(20*randn(400,3) ) ; 

octave:2> selectv = v(f ind(v( : ,2) == 7),:) 

we have selected the rows of v where the element in the second column equals 7. Now you 
see why leading computer scientists refer to octave/Matlab as "semantically suspect." It is 
a very useful language nonetheless, and you should try to learn its quirks rather than resist 
them. 

diag(v) returns the diagonal matrix with vector v as diagonal. diag(M) returns as a vector, 
the diagonal of matrix v. Thus diag(diag(v) ) is v for vector v, but diag(diag(M) ) is the 
diagonal part of matrix M. 

toeplitz(v) returns the Toeplitz matrix associated with vector v. That is 

" v(l) v(2) v(3) 
v(2) v(l) v(2) 
v(3) v(2) v(l) 



toeplitz(v) 



v(n - 1) 
v(n - 2) 



v(n) v(n-l) v(n-2) ••• v(l) 

In the more general form, toeplitz (c,r) can be used to return a assymmetric Toeplitz 
matrix. 

A matrix which is banded on the cross diagonals is evidently called a "Hankel" matrix: 

u(l) u(2) u(3) ••• u(n) 
u(2) u(3) u(4) ••• v(2) 

hankel(u,v)= u ( 3 ) u ( 4 ) u ( 5 ) v ( 3 ) 

u(n) v(2) v(3) ••• v(n) 

eig(M) returns the eigenvalues of M. [V, LAMBDA] = eig(M) returns the eigenvectors, and 
eigenvalues of M. 

kron(M,N) returns the Kronecker product of the two matrices. This is a blcok construction 
which returns a matrix where each block is an element of M as a scalar multiplied by the 
whole matrix N. 

flipud(N) flips the vector or matrix N so that its first row is last and vice versa. Similarly 
fliplr(N) flips left/right. 
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2.3 Programming and Control 

If you are going to do any serious programming in octave, you should keep your commands in a 
file, octave loads commands from '.m' files. 2 If you have the following in a file called myfunc.m: 

function [yl,y2] = myfunc(xl ,x2) 
% comments start with a "/,' 

% this function is useless, except as an example of functions. 

% input : 

% xl a number 

% x2 another number 

% output : 

% yl some output 

% y2 some output 

yl = cos(xl) . * sin(x2); 
y2 = norm(yl) ; 

then you can call this function from octave, as follows: 

octave :1> myfunc(2,3) 

ans = -0.058727 

octave :2> [a,b] = myfunc(2,3) 

a = -0.058727 

b = 0.058727 

octave :3> [a,b] = myfunc([l 2 3 4] , [1 2 3 4] ) 
a = 

0.45465 -0.37840 -0.13971 0.49468 
b = 0.78366 

Note this silly function will throw an error if xl and x2 are not of the same size. 

It is recommended that you write your functions so that they can take scalar and vector input 
where appropriate. For example, the octave builtin sine function can take a scalar and output a 
scalar, or take a vector and output a vector which is, elementwise, the sine of the input. It is not 
too difficult to write functions this way, it often only requires judicious use of . * multiplies instead 
of * multiplies. For example, if the file myfunc.m were changed to read 

yl = cos(xl) * sin(x2); 

it could easily crash if xl and x2 were vectors of the same size because matrix multiplication is not 
defined for an n x 1 matrix times another n x 1 matrix. 

An .m file does not have to contain a function, it can merely contain some octave commands. 
For example, putting the following into runner. m: 

xl = rand (4, 3) ; 

x2 = rand(size (xl) ) ; 

[a,b] = myfunc(xl ,x2) 



2 The 'm' stands for 'octave.' 
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octave allows you to call this script without arguments: 

octave :4> runner 
a = 



0.245936 
0.246414 
0.542728 
0.257607 



0.478054 
0.186454 
0.419457 
0.378558 



0.535323 
0.206279 
0.083917 
0.768188 



b = 1.3135 

octave has to know where your .m file is. It will look in the directory from which it was called. 
You can set this to something else with cd or chdir. 

You can also use the octave builtin function f eval to evaluate a function by name. For example, 
the following is a different way of calling myfunc.m: 

octave :5> [a,b] = f eval ("myfunc" ,2,3) 
a = -0.058727 
b = 0.058727 

In this form f eval seems like a way of using more keystrokes to get the same result. However, you 
can pass a variable function name as well: 

octave :6> fname = "myfunc" 
fname = myfunc 

octave :7> [a,b] = f eval (fname, 2, 3) 
a = -0.058727 
b = 0.058727 

This allows you to effectively pass functions to other functions. 



2.3.1 Logical Forks and Control 

octave has the regular assortment of 'if-then-else' and 'for' and 'while' loops. These take the 
following form: 

if exprl 

statements 
elseif expr2 

statements 
elsif expr3 

statements 

else 

statements 
end 

for var=vector 

statements 
end 

while expr 
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statements 
end 

Note that the word end is one of the most overloaded in octave/Matlab. It stands for the last 
index of a vector of matrix, as well as the exit point for for loops, if statements, switches, etc. To 
simplify debugging, it is also permissible to use endif to end an if statement, endfor to end a 
for loop, etc.. 

The test expressions may use the logical conditionals: >, <, >=, <=, ==, ~=. Do not use the 
assignment operator = as a conditional, as this can lead to disastrous results, as in C. 
Here are some examples: 

7oCompute the sign of x 
if x > 

s = 1; 
elseif x == 

s = 0; 
else 

s = -1; 
end 

%a ' regular ' for loop 
for i=l:10 

sm = sm + i; 
end 

%an 'irregular' for loop 
for i=[l 2 3 5 8 13 21 34] 

f sm = f sm + i ; 
end 

while (sin(x) > 0) 

x = x * pi; 
end 

2.4 Plotting 

Plotting is one area in which there are some noticeable differences between octave and Matlab. 
The commands and examples given herein are for octave, but the commands for Matlab are not 
too different, octave ships its plotting commands to Gnuplot. 

The main plot command is plot. You may also use semilogx, semilogy,loglog for 2D plots 
with log axes, and contour and mesh for 3D plots. Use the help command to get the specific 
syntax for each command. We present some examples: 

n = 100; 

X = pi .* ((l:n) ./ n); 
Y = sin(X) ; 
"/just plot Y 
plot(Y) ; 
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Zplot Y, but with the 'right' X axis labels 
plot(X.Y) ; 

W = sqrt(Y); 
plot (W) ; 

Zplot W, but with the 'right' X axis labels 
plot(Y.W) ; 

The output from these commands is seen in Figure 2.1. In particular, you should note the 
difference between plotting a vector, as in Figure 2.1(c) versus plotting the same vector but with 
the appropriate abscissa values, as in Figure 2.1(d). 




(c) W = W (d) W versus Y 



Figure 2.1: Four plots from octave. 

Some magic commands are required to plot to a file. For octave, I recommend the following 
magic formula, which replots the figure to a file: 

%call the plot commands before this line 
gset term postscript color; 
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gset output "filename.ps"; 
replot ; 

gset term xll; 

gset output "/dev/null" ; 

In Matlab, the commands are something like this: 

%call the plot commands before this line 
print (gcf , ' -deps ' , ' f ilename . eps ' ) ; 
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Exercises 

(2.1) What do the following pieces of octave/Matlab code accomplish? 

(a) x = (0:40) ./ 40; 

(b) a = 2; 
b = 5; 

x = a + (b-a) .* (0:40) ./ 40; 

(c) x = a + (b-a) .* (0:40) ./ 40; 
y = sin(x) ; 

plot(x,y) ; 

(2.2) Implement the naive quadratic formula to find the roots of x 2 + bx + c = 0, for real b, c. Your 
code should return 

-b ± Vfr 2 - 4c 
2 ' 

Your m-file should have header line like: 
function [xl,x2] = naivequad(b , c) 

Test your code for (b, c) = (l x 10 15 , l) . Do you get a spurious root? 

(2.3) Implement a robust quadratic formula (cf. Example Problem 1.9) to find the roots of x 2 + 
bx + c = 0. Your m-file should have header line like: 

function [xl,x2] = robustquad(b,c) 

Test your code for (b, c) = (l x 10 15 , l) . Do you get a spurious root? 

(2.4) Write octave/Matlab code to find a fixed point for the cosine, i.e., some x such that x = 
cos(x). Do this as follows: pick some initial value xq, then let Xj+i = cos(xj) for i = 
0,1,..., n. Pick n to be reasonably large, or choose some convergence criterion (i.e., ter- 
minate if |xj+i — Xi\ < 1 x 1CP 10 ). Does your code always converge? 

(2.5) Write code to implement the factorial function for integers: 
function [nfact] = factorial (n) 

where n factorial is equal to 1 ■ 2 ■ 3 ■ ■ • (n — 1) ■ n. Either use a 'for' loop, or write the function 
to recursively call itself. 



Chapter 3 

Solving Linear Systems 



A number of problems in numerical analysis can be reduced to, or approximated by, a system of 
linear equations. 



3.1 Gaussian Elimination with Naive Pivoting 

Our goal is the automatic solution of systems of linear equations: 



a n xi + 012^2 + 013X3 + 

a 2 \X\ + 022^2 + 023^3 + 

a 3i xi + 03222 + 033^3 + 



a n \x\ + a n2 X2 + a n3 x 3 + 



- (l\ n X n 

+ a2nX n 
3n X n 



6l 

b 2 
b 3 



&nnXn — b r , 



In these equations, the and bi are given real numbers. We also write this as 

Ax = b, 

where A is a matrix, whose element in the i th row and j th column is a^-, and b is a column vector, 
whose z th entry is b^. 

This gives the easier way of writing this equation: 



Oil 


ai2 


a 13 ■ ■ 


Gin 




Xi 




' h ' 


021 


022 


023 • • 


G2n 




X 2 




b 2 


031 


032 


033 • • 


• G3n 




x 3 




b 3 




«n2 


«n3 • • 


"'Till 




Xn 




b n 



3.1.1 Elementary Row Operations 

You may remember that one way to solve linear equations is by applying elementary row operations 
to a given equation of the system. For example, if we are trying to solve the given system of 
equations, they should have the same solution as the following system: 
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fit 1 


f7 1 o 




"'In 




A l 






an 


022 


023 • • 


02n 




^2 






031 




033 • • 


03n 




^3 






nan 


na i2 


na i3 ■ ■ 


K(li n 








ubi 


dnl 


a n 2 


«n3 






2-n 







where k is some given number which is not zero. It suffices to solve this system of linear 
equations, as it has the same solution(s) as our original system. Multiplying a row of the system 
by a nonzero constant is one of the elementary row operations. 

The second elementary row operation is to replace a row by the sum of that row and a constant 
times another. Thus, for example, the following system of equations has the same solution as the 
original system: 



a u 


au 


ai 3 


ain 




Xl 




' h 


an 


022 


«23 


a2n 




X2 








032 


«33 


03n 




%3 




bz 




a (i-l)2 


«(i-l)3 


a (i-l)n 




X(i-l) 




b(i-i) 


an + fiaji 


ai2 + (3aj 2 


a i3 + (3a j% ■ ■ 


ain Pajn 




x% 




h + f3b j 


a n i 


a n 2 


a n 3 


ann 








bn 



We have replaced the i row by the i row plus (3 times the j row. 
The third elementary row operation is to switch rows: 



an 


ai2 


an • • 


ain 




Xl 




' h ' 
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a n 2 


a n 3 ■ ■ 


ann 




Xn 




b n 



We have here switched the second and third rows. The purpose of this e.r.o. is mainly to make 
things look nice. 

Note that none of the e.r.o. 's change the structure of the solution vector x. For this reason, it is 
customary to drop the solution vector entirely and to write the matrix A and the vector b together 
in augmented form: 
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ai3 ■ ■ 


a In 
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The idea of Gaussian Elimination is to use the elementary row operations to put a system into 
upper triangular form then use back substitution. We'll give an example here: 
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Example Problem 3.1. Solve the set of linear equations: 

x\ + x 2 - x 3 = 2 
— 3xi — 4x 2 + 4x3 = — 7 
2xi + lx 2 + lx 3 = 7 



Solution: We start by rewriting in the augmented form: 
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-1 


2 


-3 


-4 


4 


-7 


2 


1 


1 


7 







We add 3 times the first row to the second, and —2 times the first row to the third to get: 

-1 
1 

3 

We now add —1 times the second row to the third row to get: 

-1 
1 

2 

The matrix is now in upper triangular form: there are no nonzero entries below the diagonal. This 
corresponds to the set of equations: 

x\ + x 2 - x 3 = 2 
-X2 + x 3 = -1 
2x 3 = 4 

We now solve this by back substitution. Because the matrix is in upper triangular form, we can 
solve X3 by looking only at the last equation; namely X3 = 2. However, once X3 is known, the second 
equation involves only one unknown, x 2 , and can be solved only by x 2 = 3. Then the first equation 
has only one unknown, and is solved by x\ = 1. H 

All sorts of funny things can happen when you attempt Gaussian Elimination: it may turn 
out that your system has no solution, or has a single solution (as above), or an infinite number 
of solutions. We should expect that an algorithm for automatic solution of systems of equations 
should detect these problems. 



3.1.2 Algorithm Terminology 

The method outlined above is fine for solving small systems. We should like to devise an algorithm 
for doing the same thing which can be applied to large systems of equations. The algorithm will 
take the system (in augmented form): 
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The algorithm then selects the first row as the pivot equation or pivot row, and the first element of 
the first row, an is the pivot element. The algorithm then pivots on the pivot element to get the 
system: 






V o 



an ai3 



a 



22 
J 32 



h n2 



a 23 



a 



33 



l fl3 



a\ n 


bi 


\ 


a 2n 








b' 3 





Where 



b' 




b'n ) 



(2<i<n,l<j<n) 



Effectively we are carrying out the e.r.o. of replacing the i th row by the i th row minus (f^-) ti 



times 



the first row. The quantity is the multiplier for the i th row. 

Hereafter the algorithm will not alter the first row or first column of the system. Thus, the 
algorithm could be written recursively. By pivoting on the second row, the algorithm then generates 
the system: 

/ an 





an 


Ol3 • 


a\ n 


bi 


\ 


a 22 


a 23 ' 


a 2n 


b'2 







a 33 ' 




b'l 





V o o 



l n3 



In this case 



% = 
bi = 



l 2j 



b' 



K ) 



(3 < i < n, 1 < j < n) 



3.1.3 Algorithm Problems 

The pivoting strategy we examined in this section is called 'naive' because a real algorithm is a bit 
more complicated. The algorithm we have outlined is far too rigid-it always chooses to pivot on 
the k th row during the k th step. This would be bad if the pivot element were zero; in this case all 
the multipliers — are not defined. 

Bad things can happen if a^k is merely small instead of zero. Consider the following example: 

Example 3.2. Solve the system of equations given by the augmented form: 



-0.0590 0.2372 
0.1080 -0.4348 



-0.3528 
0.6452 



Note that the exact solution of this system is x\ = 10, X2 = 1. Suppose, however, that the algorithm 
uses only 4 significant figures for its calculations. The algorithm, naively, pivots on the first 
equation. The multiplier for the second row is 

ai08 ° -1.830508..., 



-0.0590 



which will be rounded to —1.831 by the algorithm. 
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The second entry in the matrix is replaced by 

-0.4348 - (-1.831)(0.2372) = -0.4348 + 0.4343 = -0.0005, 

where the arithmetic is rounded to four significant figures each time. There is some serious sub- 
tractive cancellation going on here. We have lost three figures with this subtraction. The errors 
get worse from here. Similarly, the second vector entry becomes: 

0.6452 - (-1.831)(-0.3528) = 0.6452 - 0.6460 = -0.0008, 

where, again, intermediate steps are rounded to four significant figures, and again there is subtrac- 
tive cancelling. This puts the system in the form 



-0.0590 0.2372 
-0.0005 



-0.3528 
-0.0008 



When the algorithm attempts back substitution, it gets the value 

-0.0008 
X2 = ^00005 = 

This is a bit off from the actual value of 1. The algorithm now finds 

xi = (-0.3528 - 0.2372 • 1.6) /-0.059 = (-0.3528 - 0.3795) /-0.059 = (-0.7323) /-0.059 = 12.41, 
where each step has rounding to four significant figures. This is also a bit off. 



3.2 Pivoting Strategies for Gaussian Elimination 

Gaussian Elimination can fail when performed in the wrong order. If the algorithm selects a zero 
pivot, the multipliers are undefined, which is no good. We also saw that a pivot small in magnitude 
can cause failure. As here: 

exi +x 2 = l 
xi + x 2 = 2 

The naive algorithm solves this as 

l-X 2 1 
Xi = = 

e 1 — e 

If e is very small, then ^ is enormous compared to both 1 and 2. With poor rounding, the algorithm 
solves X2 as 1. Then it solves xi = 0. This is nearly correct for X2, but is an awful approximation 
for xi. Note that this choice of xi,X2 satisfies the first equation, but not the second. 
Now suppose the algorithm changed the order of the equations, then solved: 



Xi + X2 = 2 
exi + x 2 = 1 
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The algorithm solves this as 

1 - 2e 

X2 = — 

Xl = 2 - x 2 

There's no problem with rounding here. 

The problem is not the small entry per se: Suppose we use an e.r.o. to scale the first equation, 
then use naive G.E.: 

1 1 

Xl + -Xl = 

i € 

xi + x 2 = 2 

This is still solved as 

Xl 

and rounding is still a problem. 

3.2.1 Scaled Partial Pivoting 

The naive G.E. algorithm uses the rows 1,2,..., n-1 in order as pivot equations. As shown above, 
this can cause errors. Better is to pivot first on row £\, then row £2, etc, until finally pivoting on 
row £ n -i, for some permutation {£i}™ =1 of the integers 1, 2, . . . , n. The strategy of scaled partial 
pivoting is to compute this permutation so that G.E. works well. 

In light of our example, we want to pivot on an element which is not small compared to other 
elements in its row. So our algorithm first determines "smallness" by calculating a scale, row-wise: 

Si = max lajJ . 
i<i<« 

The scales are only computed once. 

Then the first pivot, £1, is chosen to be the i such that 

Si 

is maximized. The algorithm pivots on row £1, producing a bunch of zeros in the first column. Note 
that the algorithm should not rearrange the matrix-this takes too much work. 
The second pivot, £2, is chosen to be the i such that 

Si 

is maximized, but without choosing £2 = l\. The algorithm pivots on row £2, producing a bunch of 
zeros in the second column. 

In the k th step £k is chosen to be the i not among £i,£2,..., £k~i such that 

Si 



2-± 

e_ 

1 — X2 
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is maximized. The algorithm pivots on row £^, producing a bunch of zeros in the k column. 

The slick way to implement this is to first set £i = i for i = 1,2, ... ,n. Then rearrange this 
vector in a kind of "bubble sort": when you find the index that should be l\, swap them, i.e., find 
the j such that £j should be the first pivot and switch the values of £\,£j. 

Then at the fc th step, search only those indices in the tail of this vector: i.e., only among £j for 
k < j < n, and perform a swap. 

3.2.2 An Example 

We present an example of using scaled partial pivoting with G.E. It's hard to come up with an 
example where the numbers do not come out as ugly fractions. We'll look at a homework question. 



1 2 


-1 
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15 \ 


4 


4 
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11 
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3 
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V 6 


5 


4 


17 


31 j 



The scales are as follows: s\ = 7, s 2 = 7, S3 = 3, S4 = 17. 

We pick i\. It should be the index which maximizes |aji| /sj. These values are: 

2 4 2 6 
7' 7'3' 17' 

We pick £i=3, and pivot: 
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10 J 



We pick £2- It should not be 3, and should be the index which maximizes (a^l /si- These values 

are: 

2 2 2 
7' 7' 17' 

We have a tie. In this case we pick the second row, i.e., £2 = 2. We pivot: 
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-2 
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-3 




2 
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7 




V 





3 


7 


13 





The matrix is in permuted upper triangular form. We could proceed, but would get a zero 
multiplier, and no changes would occur. 

If we did proceed we would have £3 = 4. Then l\ = 1. Our row permutation is 3, 2,4, 1. When 
we do back substitution, we work in this order reversed on the rows, solving x±, then X3,X2,x±. 

We get X4 = 1, so 

^3 = ^ (13 -7*1) = 2 

x 2 = - (-3- 1 * 1 + 2*2) = 

xi = - (7 - 3*1-1*2-1*0) = ! 
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3.2.3 Another Example and A Real Algorithm 

Sometimes we want to solve 

Ax = b 

for a number of different vectors b. It turns out we can run G.E. on the matrix A alone and come up 
with all the multipliers, which can then be used multiple times on different vectors b. We illustrate 
with an example: 

/ 1 2 4 1 \ rr 
4 2 12 2 
2 12 3' 3 
\ 1 3 2 1 / [4 

The scale vector is s = [ 4 4 3 3] T . 

Our scale choices are |, |, |, |. We choose £\ = 2, and swap £±,£2- In the places where there 
would be zeros in the real matrix, we will put the multipliers. We will illustrate them here boxed: 



Mr 



Mi = 
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, i = 
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5 / 







Our scale choices are §,§,§• We choose £2 = 4, and so swap £2, £4: 
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Our scale choices are |g, \ . We choose ^3 = 1, and so swap ^3,^4: 



M, 



V 



27 
10 



1 \ 



5 

2 
11 

9 

1 

2 



2 
4 
1 
3 



Now suppose we had to solve the linear system for b = [ — 1 8 2 l] T - 
We scale 6 by the multipliers in order: £\ = 2, so, we sweep through the first column of M3, 
picking off the boxed numbers (your computer doesn't really have boxed variables), and scaling b 
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appropriately: 



This continues: 
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We then perform a permuted backwards substitution on the augmented system 

/ 

4 2 



This proceeds as 

X4 
£3 
^2 
Xi 



-2 9 -6 

~Tl7 ~ T7 
10/12 1 -6 

27 V T ~ 5 T7 



2 / 1-6 7 



I ' 8 



-6 

2— - x 3 - 2x 2 



Fill in your own values here. 

3.3 LU Factorization 

We examined G.E. to solve the system 
where A is a matrix: 



Ax = b, 

an «12 «13 

«21 «22 0,23 
«31 «32 ^33 



a In 
02n 
«3n 



Onl «n2 On3 ' 

We want to show that G.E. actually factors A into lower and upper triangular parts, that is A = LU, 
where 
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31 
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f"aV ^n2 ^n3 ' 1 

We call this a LU Factorization of A. 
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3.3.1 An Example 

We consider solution of the following augmented form: 
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(3.2) 



The naive G.E. reduces this to 
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We are going to run the naive G.E., and see how it is a LU Factorization. Since this is the naive 
version, we first pivot on the first row. Our multipliers are 2, 3, 1. We pivot to get 
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Careful inspection shows that we've merely multiplied A and b by a lower triangular matrix Mi: 

1 " 



-1 1 _ 

The entries in the first column are the negative e.r.o. multipliers for each row. Thus after the first 
pivot, it is like we are solving the system 

MiAz; = Mib. 



We pivot on the second row to get: 
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The multipliers are 1,-1. We can view this pivot as a multiplication by M2, with 



M 2 = 
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We are now solving 



M 2 MiAx = M 2 Mib. 



3.3. LU FACTORIZATION 



35 



We pivot on the third row, with a multiplier of —1. Thus we get 



We have multiplied by M3 : 



We are now solving 



But we have an upper triangular form, that is, if we let 
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Then we have 



M 3 M 2 MiA= U, 



A= (M 3 M 2 Mi) _1 U, 



A 
A 



Mi" 
LU. 



• 1 M 2 " 1 M 3 - 1 U, 



We are hoping that L is indeed lower triangular, and has ones on the diagonal. It turns out that 
the inverse of each Mj matrix has a nice form (See Exercise (3.6)). We write them here: 
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This is really crazy: the matrix L looks to be composed of ones on the diagonal and multipliers 
under the diagonal. 

Now we check to see if we made any mistakes: 
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3.3.2 Using LU Factorizations 

We see that the G.E. algorithm can be used to actually calculate the LU factorization. We will look 
at this in more detail in another example. We now examine how we can use the LU factorization 
to solve the equation 

Ax = b, 

Since we have A = LU, we first solve 

Iz = b, 

then solve 

\Jx = z. 

Since L is lower triangular, we can solve for z with a forward substitution. Similarly, since U is 
upper triangular, we can solve for x with a back substitution. We drag out the previous example 
(which we never got around to solving): 
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We had found the LU factorization of A as 



A = 



So we solve 



We get 



Now we solve 



We get the ugly solution 
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3.3.3 Some Theory 

We aren't doing much proving here. The following theorem has an ugly proof in the Cheney & 
Kincaid [7]. 

Theorem 3.3. If A is an n x n matrix, and naive Gaussian Elimination does not encounter a zero 
pivot, then the algorithm generates a LU factorization of A, where L is the lower triangular part of 
the output matrix, and U is the upper triangular part. 

This theorem relies on us using the fancy version of G.E., which saves the multipliers in the 
spots where there should be zeros. If correctly implemented, then, L is the lower triangular part 
but with ones put on the diagonal. 

This theorem is proved in Cheney & Kincaid [7]. This appears to me to be a case of something 
which can be better illustrated with an example or two and some informal investigation. The proof 
is an unillustrating index-chase-read it at your own risk. 

3.3.4 Computing Inverses 

We consider finding the inverse of A. Since 

AA" 1 = I, 

then the j th column of the inverse A -1 solves the equation 

Ace — c j . 

where ej is the column matrix of all zeros, but with a one in the j th position. 

Thus we can find the inverse of A by running n linear solves. Obviously we are only going 
to run G.E. once, to put the matrix in LU form, then run n solves using forward and backward 
substitutions. 

3.4 Iterative Solutions 

Recall we are trying to solve 

Ax = b. 

We examine the computational cost of Gaussian Elimination to motivate the search for an alter- 
native algorithm. 

3.4.1 An Operation Count for Gaussian Elimination 

We consider the number of floating point operations ("flops") required to solve the system Ate = b. 
Gaussian Elimnation first uses row operations to transform the problem into an equivalent problem 
of the form \Jx = b', where U is upper triangular. Then back substitution is used to solve for x. 
First we look at how many floating point operations are required to reduce 
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First a multiplier is computed for each row. Then in each row the algorithm performs n 
multiplies and n adds. This gives a total of (n— 1) + (n — l)n multiplies (counting in the computing 
of the multiplier in each of the (n — 1) rows) and (n — l)n adds. In total this is 2n 2 — n — 1 floating 
point operations to do a single pivot on the n by n system. 

Then this has to be done recursively on the lower right subsystem, which is an (n— 1) by (n — 1) 
system. This requires 2(n — 1) 2 — (n — 1) — 1 operations. Then this has to be done on the next 
subsystem, requiring 2(n — 2) 2 — (n — 2) — 1 operations, and so on. 

In total, then, we use /„ total floating point operations, with 

n n n 

j=l j=l j=l 



Recalling that 



n ^ n ^ 

^J 2 = ^W(n + l)(2n + l), and ^j = -(n)(n + l), 



j'=i 



We find that 



1 2 o 

/„ = -(4n - l)n(n + 1) - n w -n 6 . 
b o 



Now consider the costs of back substitution. To solve 
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for x n requires only a single division. Then to solve for x n -\ we compute 

1 



Xn—l 



Q-n—l.n—l 



[bn— 1 Q"n— l,n%n] > 



and requires 3 flops. Similarly, solving for x n _2 requires 5 flops. Thus in total back substitution 
requires B n total floating point operations with 



B n = 2j — 1 = n(n — 1) — n = n(n — 2)mn z 

3=1 
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3.4.2 Dividing by Multiplying 

We saw that Gaussian Elimination requires around |ra 3 operations just to find the LU factorization, 
then about n 2 operations to solve the system, when A is n x n. When n is large, this may take too 
long to be practical. Additionally, if A is sparse (has few nonzero elements per row), we would like 
the complexity of our computations to scale with the sparsity of A. Thus we look for an alternative 
algorithm. 

First we consider the simplest case, n = 1. Suppose we are to solve the equation 



Ax = b. 



for scalars A, b. We solve this by 



\ 1 1 

A ~ ujA ~ 1 - (1 - ojA) 



uib = 



1-r 



uib, 



where uo / is some real number chosen to "weight" the problem appropriately, and r = 1 — to A. 
Now suppose that oj is chosen such that \r\ < 1. This can be done so long as A / 0, which would 
have been a problem anyway. Now use the geometric expansion: 



1 + r + r 2 + r 6 + . . . 



Because of the assumption \r\ < 1, the terms r n converge to zero as n 
approximate solution to our one dimensional problem as 



oo. This gives the 



1 + r + r 2 + r 3 + . . . + r k 



Ldb + 

Lob + r 



r + r 2 +r 3 + ...+r k 



ujb 
ub 



1 + r + r 2 + 



+ r 



k-l 



cob 



This suggests an iterative approach to solving Ax = b. First let x^ = cob, then let 

x ( fc ) = ub + rx (fc_1) . 

The iterates i^' will converge to the solution of Ax = b if \r\ < 1. 

You should now convince yourself that because r n — > 0, that the choice of the initial iterate x^ ^ 
was immaterial, i.e., that under any choice of initial iterate convergence is guaranteed. 

We now translate this scalar result into the vector case. The algorithm proceeds as follows: first 
fix some initial estimate of the solution, x^°\ A good choice might be cob, but this is not necessary. 
Then calculate successive approximations to the actual solution by updates of the form 



x 



W =U jb+ (I -ujk)x 



(k-l) 



It turns out that we can consider a slightly more general form of the algorithm, one in which 
successive iterates are defined implicitly. That is we consider iterates of the form 



Qa;( fc+1 ) = (Q - ljA) ajW + cob, 



(3.3) 



for some matrix Q, and some scaling factor to. Note that this update relies on vector additions 
and possibly by premultiplication of a vector by A or Q. In the case where these two matrices are 
sparse, such an update can be relatively cheap. 
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Now suppose that as k — > oo, x^ converges to some vector x* , which is a fixed point of the 
iteration. Then 

Qx* = (Q -ujA)x* + ujb, 
Qx* = Qx* - ujAx* + ujb, 
u)Ax* = oob, 
Ax* = b. 

We have some freedom in choosing Q, but there are two considerations we should keep in mind: 

1. Choice of Q affects convergence and speed of convergence of the method. In particular, we 
want Q to be similar to A. 

2. Choice of Q affects ease of computing the update. That is, given 

z = (Q - A) + b, 
we should pick Q such that the equation 

Qa;(fc+i) = z 

is easy to solve exactly. 

These two goals conflict with each other. At one end of the spectrum is the so-called "impossible 
iteration," at the other is the Richardsons. 

3.4.3 Impossible Iteration 

I made up the term "impossible iteration." But consider the method which takes Q to be A. This 
seems to be the best choice for satisfying the first goal. Letting uj = 1, our method becomes 

Aa ,(fc+i) = (A — A) £c( fe ) + b = b. 

This method should clearly converge in one step. However, the second goal is totally ignored. 
Indeed, we are considering iterative methods because we cannot easily solve this linear equation in 
the first place. 

3.4.4 Richardson Iteration 

At the other end of the spectrum is the Richardson Iteration, which chooses Q to be the identity 
matrix. Solving the system 

Q^ +1 ) = z 

is trivial: we just have x^ k+1 ^ = z. 

Example Problem 3.4. Use Richardson Iteration with w = lon the system 





" 6 


1 


1 " 




" 12 " 


A = 


2 


4 





b = 







1 


2 


6 




6 



Solution: We let 





" l 





" 




" -5 


-1 


-1 


Q = 





1 





, (Q-A) = 


-2 


-3 













1 




-1 


-2 


-5 
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We start with an arbitrary x^\ say = [2 2 2] T . We get x^ = [-2 - 10 - 10] T , and x^ = 
[42 34 78] T . 

Note the real solution is x = [2 — 1 1] T . The Richardson Iteration does not appear to converge 
for this example, unfortunately. H 

Example Problem 3.5. Apply Richardson Iteration with w = 1/6 on the previous system. 
Solution: Our iteration becomes 








-1/6 


-1/6 " 




' 2 " 


x (k+l) = 


-1/3 


1/3 





x^ + 







_ -1/6 


-1/3 







1 



We start with the same x^ as previously, x^ = [2 2 2] T . We get x^ = [4/3 0] T , x^ = 
[2 - 4/9 7/9] T , and finally x^ = [2 - 0.99998 0.99998] T . 

Thus, the choice of u has some affect on convergence. H 

We can rethink the Richardson Iteration as 

X ( k +V = (| _ W A) cc^ + cub = + w (b - AajW) . 

Thus at each step we are adding some scaled version of the residual, defined as b — Ax , to the 
iterate. 



3.4.5 Jacobi Iteration 



The Jacobi Iteration chooses Q to be the matrix consisting of the diagonal of A. This is more 
similar to A than the identity matrix, but nearly as simple to invert. 

Example Problem 3.6. Use Jacobi Iteration, with u> = 1, to solve the system 
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1 " 








" 12 " 














A = 
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4 
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b 





















1 
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6 








6 






We let 




























" 6 





" 













-1 


-1 " 




- 1 

6 





Q = 





4 





, (Q- 


-A) 






2 










Q- l = 





1 

4 










6 










1 


-2 














We start with an arbitrary say x^ = [2 2 2] 

[t ~ I T ] T • Continuing, we find that x^ w [1.987 - 1.019 0.981] 
Note the real solution is x = [2 — 1 1] T . 



We get x^ = [§ - 1 0] T . Then x^ = 

T 



There is an alternative way to describe the Jacobi Iteration for uj = 1. By considering the update 
elementwise, we see that the operation can be described by 



x 



(fe+i) _ 1 




Thus an update takes less than 2n 2 operations. In fact, if A is sparse, with less than k nonzero 
entries per row, the update should take less than 2nk operations. 
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3.4.6 Gauss Seidel Iteration 

The Gauss Seidel Iteration chooses Q to be lower triangular part of A, including the diagonal. In 
this case solving the system 

Qx^ = z 

is performed by forward substitution. Here the Q is more like A than for Jacobi Iteration, but 
involves more work for inverting. 

Example Problem 3.7. Use Gauss Seidel Iteration to again solve for 





" 6 
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1 " 




" 12 " 


A = 
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4 





b = 
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2 


6 




6 



Solution: We let 



Q 



6 
2 4 
1 2 6 



(Q-A) 



-1 





-1 






We start with an arbitrary x(°), say x<W = [2 2 2]. We get xW = [§ - § l] . Then x 

[35 _ 35 il T 
L 18 36 1 1 ■ 

Already this is fairly close to the actual solution x = [2 — 1 1] T . 



,(2) 



Just as with Jacobi Iteration, there is an easier way to describe the Gauss Seidel Iteration. In 
this case we will keep a single vector x and overwrite it, element by element. Thus for j = 1, 2, . . . , n, 
we set 

Xj *~ I °j 2s a i* x * I • 

" \ i=l&j j 

This looks exactly like the Jacobi update. However, in the sum on the right there are some "old" 
values of Xj and some "new" values; the new values are those X{ for which i < j. 

Again this takes less than 2n 2 operations. Or less than Ink if A is sufficiently sparse. 

An alteration of the Gauss Seidel Iteration is to make successive "sweeps" of this redefinition, 
one for j = 1, 2, . . . , n, the next for j = n, n — 1, . . . , 2, 1. This amounts to running Gauss Seidel 
once with Q the lower triangular part of A, then running it with Q the upper triangular part. This 
iterative method is known as "red-black Gauss Seidel." 

3.4.7 Error Analysis 

Suppose that x is the solution to equation 3.4. Define the error vector: 
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Now notice that 



x 



x 



X 



(fc+1) = 
(fc+1) = 

(fc+1) _ 



Q _1 (Q - ujA) £c (fe) + Q _1 u;b, 
Q-^ajW - wQ _1 AajW + wQ^Aa;, 



a? 



(fc+i) 



a; = a; 



(fc) 



IE 



(fc) 



a; 



- x - wQ^A (x^ - x^j , 

(fc) 



= (l-wQ^A) cW. 



e (fe+i) 

Reusing this relation we find that 

e« = (I - wQ^A) e^" 1 ), 

= (l-cQ- 1 A) fc e(°). 

We want to ensure that e^ k+1 ^ is "smaller" than e^ k \ To do this we recall matrix and vector norms 
from Subsection 1.4.1. 



,(*) 



(l-wQ- 1 A) fc c(°) < j|l -t^Q^A 



,(o) 



(See Example Problem 1.29.) 

Thus our iteration converges (e^ goes to the zero vector, i.e., x^ — ► x) if 



|l -wQ _1 A|| 2 < 1. 



This gives the theorem: 



Theorem 3.8. An iterative solution scheme converges for any starting x^ if and only if all 
eigenvalues of I — wQ _1 A are less than 1 in absolute value, i.e., if and only if 

||l -<jQ _1 A|| 2 < 1 

Another way of saying this is "the spectral radius of I — wQ _1 A is less than 1." 

In fact, the speed of convergence is decided by the spectral radius of the matrix-convergence 
is faster for smaller values. Recall our introduction to iterative methods in the scalar case, where 
the result relied on uj being chosen such that |1 — ujA\ < 1. You should now think about how 
eigenvalues generalize the absolute value of a scalar, and how this relates to the norm of matrices. 

Let y be an eigenvector for Q _1 A, with corresponding eigenvalue A. Then 

(I - wQ _1 A) y = y- ujQ- l Ay = y-uXy = (l- uX) y. 

This relation may allow us to pick the optimal uj for given A,Q. It can also show us that 
sometimes no choice of uj will give convergence of the method. There are a number of different 
related results that show when various methods will work for certain choices of uj. We leave these 
to the exercises. 
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Example Problem 3.9. Find conditions on uj which guarantee convergence of Richardson's Iter- 
ation for finding approximate iterative solutions to the system Ax = b, where 





" 6 


1 


1 " 




" 12 " 


A = 


2 


4 





b = 







1 


2 


6 




6 



Solution: By Theorem 3.8, with Q the identity matrix, we have convergence if and only if 

||l - wA|| 2 < 1 

We now use the fact that "eigenvalues commute with polynomials;" that is if f(x) is a polynomial 
and A is an eigenvalue of a matrix A, then /(A) is an eigenvalue of the matrix /(A). In this case the 
polynomial we consider is f(x) = x°—lux 1 . Using octave or Matlab you will find that the eigenvalues 
of A are approximately 7.7321,4.2679, and 4. Thus the eigenvalues of I — uiA are approximately 



1 - 7.7321w, 1 - 4.2679cj, 1 - 4u. 



With some work it can be shown that all three of these values will be less than one in absolute 
value if and only if 

< uj < — ^ — « 0.388 
7.7321 

(See also Exercise (3.10).) 

Compare this to the results of Example Problem 3.4, where for this system, ui = 1 apparently did 
not lead to convergence, while for Example Problem 3.5, with uj = 1/6, convergence was observed. 

H 



3.4.8 A Free Lunch? 

The analysis leading to Theorem 3.8 leads to an interesting possible variant of the iterative scheme. 
For simplicity we will only consider an alteration of Richardson's Iteration. In the altered algorithm 
we presuppose the existence, via some oracle, of a sequence of weightings, u>i, which we use in each 
iterative update. Thus our algorithm becomes: 

1. Select some initial iterate x^°\ 

2. Given iterate £c( fc_1 ), define 

X W = (\ - uj k A) x^ + u; k b. 
Following the analysis for Theorem 3.8, it can be shown that 

e ( fc ) = (I - u k A) e^ -1 ) 

where, again, = x^ — x, with x the actual solution to the linear system. Expanding e^ -1 ) 
similarly gives 

e W = ( |_^ A )(l^_ 1 A)e^ 2 ) 
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Remember that we want to be small in magnitude, or, better yet, to be zero. One way to 
guarantee that is zero, regardless of the choice of e(°) it to somehow ensure that the matrix 

k 

B = H I - u>iA 
1=1 

has all zero eigenvalues. 

We again make use of the fact that eigenvalues "commute" with polynomials to claim that if 
Xj is an eigenvalue of A, then 

k 

n 1 - ^ 

i=i 

is an eigenvalue of B. This eigenvalue is zero if one of the uj{ for 1 < % < k is 1/Aj. This suggests 
how we are to pick the weightings: let them be the inverses of the eigenvalues of A. 

In fact, if A has a small number of distinct eigenvalues, say m eigenvalues, then convergence to 
the exact solution could be guaranteed after only m iterations, regardless of the size of the matrix. 

As you may have guessed from the title of this subsection, this is not exactly a practical 
algorithm. The problem is that it is not simple to find, for a given arbitrary matrix A, one, some, 
or all its eigenvalues. This problem is of sufficient complexity to outweigh any savings to be gotten 
from our "free lunch" algorithm. 

However, in some limited situations this algorithm might be practical if the eigenvalues of A 
are known a priori. 
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Exercises 



(3.1) Find the LU decomposition of the following matrices, using naive Gaussian Elimination 



(a) 



3 -1 -2 
9 -1 -4 
-6 10 13 



(b) 



24 
12 
13 



16 
11 
19 



-3 
1 

-4 



6 
-2 
5 -8 



(3.2) Perform back substitution to solve the equation 



1 


3 


5 


3 " 




' -1 " 





2 


4 


3 




-1 










x = 










2 


1 


1 











2 




2 



(3.3) Perform Naive Gaussian Elimination to prove Cramer's rule for the 2D case. That is, prove 
that the solution to 



a b 




X 




" / " 


c d 




. v _ 




. a _ 



is given by 



det 



det 



a f 

c g 

a b 

c d 



and x 



det 



det 



/ b 

9 d 

a b 

c d 



(3.4) 



(3.5) 



(3.6) 



Implement Cramer's rule to solve a pair of linear equations in 2 variables. Your m-file should 
have header line like: 
function x = cramer2(A,b) 

where A is a 2 x 2 matrix, and b and x are 2x1 vectors. Your code should find the x such 
that Ax = b. (See Exercise (3.3)) 

Test your code on the following (augmented) systems: 



(a) 
(d) 



3 -2 

4 -3 
-0.0590 

0.1080 



(b) 



0.2372 
-0.4348 



1.24 
-0.744 
-0.3528 
0.6452 



-3.48 
2.088 



(c) 



1.24 -3.48 
-0.744 2.088 



1 

-0.6 



Given two lines parametrized by f(t) = at + b, and g(s) = cs + d, set up a linear 2x2 system 
of equations to find the t, s at the point of intersection of the two lines. If you were going 
to write a program to detect the intersection of two lines, how would you detect whether 
they are parallel? How is this related to the form of the solution to a system of two linear 
equations? (See Exercise (3.3)) 
Prove that the inverse of the matrix 



1 ••• 

02 1 ■■■ 

3 1 •■■ 

a n ■■■ 1 
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is the matrix 



1 

-a 2 1 

-a 3 1 

-Or, 







1 



(Hint: Multiply them together.) 
(3.7) Under the strategy of scaled partial pivoting, which row of the following matrix will be the 
first pivot row? 



10 


17 


-10 


0.1 


0.9 


-3 


3 


-3 


0.3 


-4 


0.3 


0.1 


0.01 


-1 


0.5 


2 


3 


4 


-3 


5 


10 


100 


1 


0.1 






(3.8) Let A be a symmetric positive definite n x n matrix with n distinct eigenvalues. Letting 



y(°) = 5/ ||6|| 2 5 consider the iteration 



y 



(fc+i) _ 



WW 



(a) What is ||y (fc) || 2 ? 

(b) Show that y( fc ) 

(c) Show that as k 



A k b/\\A k b\\ 2 . 

oo, y^> converges to the (normalized) eigenvector associated with the 



largest eigenvalue of A. 
(3.9) Consider the equation 



Letting x^ = [11 0] T , find the iterate a^ 1 ) by one step of Richardson's Method. And by 
one step of Jacobi Iteration. And by Gauss Seidel. 
(3.10) Let A be a symmetric n x n matrix with eigenvalues in the interval [a, with < a < (3, 
and a + (3 ^ 0. Consider Richardson's Iteration 



1 


3 


5 




' -5 " 


-2 


2 


4 


x = 


-6 


4 


-3 


-4 




10 



X 



(fc+1) 



(I — uA) x^> + ub. 



Recall that e( fc+1 ) = (I - ujA) . 

(a) Show that the eigenvalues of I — ujA are in the interval [1 — u>(3, 1 — wa]. 

(b) Prove that 

max {|A| : 1 — tu(3 < A < 1 — iva} 



ujI3 



is minimized when we choose ui such that 1 
look at the graph of something versus to.) 

(c) Show that this relationship is satisfied by u> = 2/ (a + 0). 

(d) For this choice of u show that the spectral radius of I — ooA is 

\a + [3\' 



(1 — ua) . (Hint: It may help to 
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(3.11 



(e) Show that when < a, this quantity is always smaller than 1. 

(f) Prove that if A is positive definite, then there is an uj such that Richardson's Iteration 
with this lo will converge for any choice of x (°) . 

(g) For which matrix do you expect faster convergence of Richardson's Iteration: Ai with 
eigenvalues in [10,20] or A2 with eigenvalues in [1010, 1020]? Why? 

I Implement Richardson's Iteration to solve the system Ax = b. Your m-file should have 
header line like: 

function xk = richardsons(A,b,xO,w,k) 

Your code should return x^ based on the iteration 



x 



O'+i) 



X 



(i) 



UJ 



(Ax® 



Let w take the place of uj, and let xO be the initial iterate x^°\ Test your code for A, b for 
which you know the actual solution to the problem. (Hint: Start with A and the solution x 
and generate b.) Test your code on the following matrices: 

• Let A be the Hilbert Matrix. This is generated by the octave command A = hilb(10), 
or whatever (smallish) integer. Try different values of uj, including uj = 1. 

• Let A be a Toeplitz matrix of the form: 



A 



These can be generated by the octave command A = toeplitz ([-2 100000 0]). 

Try different values of uj, including uj = —1/2. 
(3.12) Let A be a nonsingular n x n matrix. We wish to solve Ax = b. Let x^ be some starting 
vector, let Vk be span{r^,Ar^, . . . , A fc r^} , and let Vk be the set of polynomials, p(x) of 
degree k with p(0) = 1. 

Consider the following iterative method: Let cc( fe+1 ) be the x that solves 



-2 


1 


•• 


• 


1 


-2 


1 •• 


• 





1 


-2 •• 


• 








•• 


• -2 



mm 



\b - Ax\ 



2 • 



Let rW = b- Ax^ k \ 

(a) Show that if x € a;( ) + V k , then b - Ax = p(A)r^ for some p G V k . 

(b) Prove that, conversely, for any p G Vk there is some x G x(°) + Vk, such that b 
p(A)r(°). 

(c) Argue that 



Ax 



,(fc+i) 



mm 

pev k 



p(A)r 



(0) 



(d) Prove that this iteration converges in at most n steps. (Hint: Argue for the existence 
of a polynomial in V n that vanishes at all the eigenvalues of A. Use this polynomial to 
show that ||r( n )|L < 0.) 



Chapter 4 

Finding Roots 



4.1 Bisection 

We are now concerned with the problem of finding a zero for the function f(x), i.e., some c such 
that /(c) = 0. 

The simplest method is that of bisection. The following theorem, from calculus class, insures 
the success of the method 

Theorem 4.1 (Intermediate Value Theorem). If f(x) is continuous on [a,b] then for any y 
such that y is between f(a) and f(b) there is a c G [a, b] such that f(c) = y. 

The IVT is best illustrated graphically. Note that continuity is really a requirement here-a 
single point of discontinuity could ruin your whole day, as the following example illustrates. 

Example 4.2. The function f(x) = \ is not continuous at 0. Thus if G [a, b] , we cannot apply 
the IVT. In particular, if G [a,b] it happens to be the case that for every y between f(a),f(b) 
there is no c G [a, b] such that /(c) = y. 

In particular, the IVT tells us that if f(x) is continuous and we know a, b such that /(a), f(b) 
have different sign, then there is some root in [a,b]. A decent estimate of the root is c = 
We can check whether /(c) = 0. If this does not hold then one and only one of the two following 
options holds: 

1. /(a), /(c) have different signs. 

2. /(c), f(b) have different signs. 

We now choose to recursively apply bisection to either [a,c] or [c, b], respectively, depending on 
which of these two options hold. 

4.1.1 Modifications 

Unfortunately, it is impossible for a computer to test whether a given black box function is contin- 
uous. Thus malicious or incompetent users could cause a naively implemented bisection algorithm 
to fail. There are a number of easily conceivable problems: 

1. The user might give f,a,b such that f(a),f(b) have the same sign. In this case the function 
/ might be legitimately continuous, and might have a root in the interval [a, b] . If, taking 
c = /(a), f(b), /(c) all have the same sign, the algorithm would be at an impasse. We 
should perform a "sanity check" on the input to make sure f(a),f(b) have different signs. 
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2. The user might give f,a,b such that / is not continuous on [a, b], moreover has no root in 
the interval [a, b] . For a poorly implemented algorithm, this might lead to an infinite search 
on smaller and smaller intervals about some discontinuity of /. In fact, the algorithm might 
descend to intervals as small as machine precision, in which case the midpoint of the interval 
will, due to rounding, be the same as one of the endpoints, resulting in an infinite recursion. 

3. The user might give / such that / has no root c that is representable in the computer's 
memory. Recall that we think of computers as storing numbers in the form ±r x 10 fc ; given 
a finite number of bits to represent a number, only a finite number of such numbers can be 
represented. It may legitimately be the case that none of them is a root to /. In this case, 
the behaviour of the algorithm may be like that of the previous case. A well implemented 
version of bisection should check the length of its input interval, and give up if the length is 
too small, compared to machine precision. 

Another common error occurs in the testing of the signs of /(a), f(b). A slick programmer might 
try to implement this test in the pseudocode: 

if (f (a)f (b) > 0) then . . . 

Note however, that |/(a)| , |/(&)| might be very small, and that f(a)f(b) might be too small to 
be representable in the computer; this calculation would be rounded to zero, and unpredictable 
behaviour would ensue. A wiser choice is 

if (sign(f(a)) * sign(f(b)) > 0) then ... 

where the function sign (x) returns —1,0, 1 depending on whether x is negative, zero, or positive, 
respectively. 

The pseudocode bisection algorithm is given as Algorithm 1. 
4.1.2 Convergence 

We can see that each time recursive_bisection (/, a,b, . . .) is called that \b — a\ is half the length 
of the interval in the previous call. Formally call the first interval [ao,&oL an d the first midpoint 
Co- Let the second interval be [ai,&i], etc. Note that one of a±,bi will be cp, and the other will be 
either ao or bo. We are claiming that 

, _ bn-i — a n -i 

bp - ap 
~ 2 n 

Theorem 4.3 (Bisection Method Theorem). If f(x) is a continuous function on [a, b] such that 
f(a)f(b) < 0, then after n steps, the algorithm run_bisection will return c such that \c — c'\ < 
^ft^, where d is some root of /. 

4.2 Newton's Method 

Newton's method is an iterative method for root finding. That is, starting from some guess at the 
root, xp, one iteration of the algorithm produces a number x\, which is supposed to be closer to a 
root; guesses X2,x%, . . . ,x n follow identically. 

Newton's method uses "linearization" to find an approximate root. Recalling Taylor's Theorem, 
we know that 

f(x + h)^f(x) + f'(x)h. 
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Algorithm 1: Algorithm for finding root by bisection. 
Input: a function, two endpoints, a x-tolerance, and a y-tolerance 
Output: a c such that \f(c)\ is smaller than the y-tolerance. 
run_bisection(/, a, b, 5, e) 

(1) Let /a <- sign (/(a)). 

(2) Let fb<- sign (/(&)). 

(3) if /a/6 > 

(4) throw an error. 

(5) else if fafb = 

(6) if fa = 

(7) return a 

(8) else 

(9) return b 

(10) Let c ^ ^ 

(11) while 6-a> 

(12) Let /c <- /(c). 

(13) if |/c| < e 

(14) return c 

(15) if sign (fa) sign (/c) < 

(16) Let b^c, /&<-/c, c ^- s±£. 

(17) else 

(18) Let a <- c, fa <- /c, c «- 

(19) return c 



This approximation is better when /"(•) is "well-behaved" between x and x + h. Newton's method 
attempts to find some h such that 



This is easily solved as 



= f(x + h) = f(x) + f'(x)h. 
-fix) 



h 



fix) 



An iteration of Newton's method, then, takes some guess Xk and returns Xk+i defined by 



Xk+l = x k 



f(Xk) 

f'(x k y 



(4.1) 



An iteration of Newton's method is shown in Figure 4.1, along with the linearization of f(x) at 



Xk- 



4.2.1 Implementation 

Use of Newton's method requires that the function fix) be differentiable. Moreover, the derivative 
of the function must be known. This may preclude Newton's method from being used when f(x) 
is a black box. As is the case for the bisection method, our algorithm cannot explicitly check for 
continuity of f(x). Moreover, the success of Newton's method is dependent on the initial guess 
xq. This was also the case with bisection, but for bisection there was an easy test of the initial 
interval-i.e., test if /(ao)/(&o) < 0- 
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Figure 4.1: One iteration of Newton's method is shown for a quadratic function f(x). The lin- 
earization of f{x) at Xk is shown. It is clear that xt+i is a root of the linearization. It happens to 
be the case that \f(xk+i)\ is smaller than , i.e., x^+i is a better guess than xp.. 



Our algorithm will test for goodness of the estimate by looking at • The algorithm will 

also test for near-zero derivative. Note that if it were the case that f'( x k) = then h would be ill 
defined. 

Algorithm 2: Algorithm for finding root by Newton's Method. 

Input: a function, its derivative, an initial guess, an iteration limit, and a tolerance 

Output: a point for which the function has small value. 

run_newton(/, /', xO, N, tol) 



(1) 


Let x <— xO, n <— 0. 


(2) 


while n < N 


(3) 


Let fx <- f{x). 


(4) 


if \fx\ < tol 


(5) 


return x. 


(6) 


Let fpx <— f'(x). 


(7) 


if \fpx\ < tol 


(8) 


Warn u f'(x) is small 


(9) 


return x. 


(10) 


Let x <— x — fx/ fpx. 


(11) 


Let n <— n + 1. 
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4.2.2 Problems 

As mentioned above, convergence is dependent on /(x), and the initial estimate xq. A number of 
conceivable problems might come up. We illustrate them here. 

Example 4.4. Consider Newton's method applied to the function /(x) = x-? with j > 1, and with 
initial estimate Xq ^ 0. 

Note that f(x) = has the single root x = 0. Now note that 

4 ^ 1 



Xk+l =X k 73T = 1 - - )x k . 

jx J k V 3 J 

Since the equation has the single root x = 0, we find that x& is converging to the root. However, 
it is converging at a rate slower than we expect from Newton's method: at each step we have a 
constant decrease of 1 — (l/j) , which is a larger number (and thus worse decrease) when j is larger. 

Example 4.5. Consider Newton's method applied to the function f(x) = with initial estimate 
x = 3. 

Note that /(x) is continuous on R + . It has a single root at x = 1. Our initial guess is not too far 
from this root. However, consider the derivative: 

,,, . x\ - lnx 1 - lnx 

fix) = ~^— 2 = — 

x z x z 

If x > e 1 , then 1 — lnx < 0, and so /'(x) < 0. However, for x > 1, we know /(x) > 0. Thus taking 

Xk+l =X k - -J7-, r > X k . 

f'(Xk) 

The estimates will "run away" from the root x = 1. 

Example 4.6. Consider Newton's method applied to the function f(x) = sin (x) for the initial 
estimate xo 7^ 0, where xq has the odious property 2xo = tanxo- 

You should verify that there are an infinite number of such xq. Consider the identity of X\ : 

f(x ) sin(x ) 

Xl = X - — = Xq = Xq- tanx = X - 2x = -X . 

f'(x ) cos(xo) 

Now consider X2 : 

/(xi) sin(— xo) sin(xo) 

x 2 = xi - — — - = -x 7 r - x H 7 — r = -x + tan x = -x + 2x = £o- 

/'(xi) cos(-xo) cos(x ) 

Thus Newton's method "cycles" between the two values xo, — xo- 

Of course, Newton's method may find some iterate Xk for which f'(x k ) = 0, in which case, there 
is no well-defined x^+i. 



4.2.3 Convergence 

When Newton's Method converges, it actually displays quadratic convergence. That is, if e k = x^— r, 

1 1 2 

where r is the root that the Xk are converging to, that |efc+i| < C | | . If, for example, it were 
the case that C = 1, then we would double the accuracy of our root estimate with each iterate. 
That is, if eo were 0.001, we would expect e\ to be on the order of 0.000001. The following theorem 
formalizes our claim: 
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Theorem 4.7 (Newton's Method Convergence). If f(x) has two continuous derivatives, and 
r is a simple root of f(x), then there is some D such that if \xq — r\ < D, Newton's method will 
converge quadratically to r. 

The proof is based on arguments from real analysis, and is omitted; see Cheney &; Kincaid for 
the proof [7]. Take note of the following, though: 

1. The proof requires that r be a simple root, that is that f'(r) ^ 0. When this does not hold 
we may get only linear convergence, as in Example 4.4. 

2. The key to the proof is using Taylor's theorem, and the definition of Newton's method to 
show that 

k+1 2f'(x k ) ' 

where is between Xk and r = Xk + e&. 

The proof then proceeds by showing that there is some region about r such that 

(a) in this region \ f"{x)\ is not too large and is not too small, and 

(b) if Xk is in the region, then Xk+i is in the region. 

In particular the region is chosen to be so small that if Xk is in the region, then the factor 
e\ will outweigh the factor /^\f'( x k)\ ■ You can roughly think of this region as an 

"attractor" region for the root. 

3. The theorem never guarantees that some Xk will fall into the "attractor" region of a root r of 
the function, as in Example 4.5 and Example 4.6. The theorem that follows gives sufficient 
conditions for convergence to a particular root. 

Theorem 4.8 (Newton's Method Convergence II [2]). If f(x) has two continuous derivatives 
on [a,b], and 

1. f(a)f(b) < 0, 

2. /'(x)^0on [a, b], 

3. fix) does not change sign on [a, b], 

4. Both |/(a)| < (6 - a) |/'(a)| and |/(6)| < (b - a) \f'(b)\ hold, 

then Newton's method converges to the unique root of f(x) = for any choice of xq £ [a, b] . 

4.2.4 Using Newton's Method 

Newton's Method can be used to program more complex functions using only simple functions. 
Suppose you had a computer which could perform addition, subtraction, multiplication, division, 
and storing and retrieving numbers, and it was your task to write a subroutine to compute some 
complex function g(-). One way of solving this problem is to have your subroutine use Newton's 
Method to solve some equation equivalent to g(z) — x = 0, where z is the input to the subroutine. 
Note that it is assumed the subroutine cannot evaluate g(z) directly, so this equation needs to be 
modified to satisfy the computer's constraints. 

Quite often when dealing with problems of this type, students make the mistake of using New- 
ton's Method to try to solve a linear equation. This should be an indication that a mistake was 
made, since Newton's Method can solve a linear equation in a single step: 

Example 4.9. Consider Newton's method applied to the function f(x) = ax + b. The iteration is 
given as 

axk + b 

x k+ i <- x k • 

a 

This can be rewritten simply as Xk+i < b/a. 
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The following example problems should illustrate this process of "bootstrapping" via Newton's 
Method. 

Example Problem 4.10. Devise a subroutine using only subtraction and multiplication that can 
find the multiplicative inverse of an input number z, i.e., can find (1/z) . 

Solution: We are tempted to use the linear function f(x) = (1/z) — x. But this is a linear equation 
for which Newton's Method would reduce to Xk+i (1/z) • Since the subroutine can only use 
subtraction and multiplication, this will not work. 

Instead apply Newton's Method to f(x) = z — (1/x) . The Newton step is 

z — (1/xk) 2 

x k+l ^~ x k 7Z i 2\ = x k ~ ZX k + Xk = x k (2 — ZXk) ■ 

Note this step uses only multiplication and subtraction. The subroutine is given in Algorithm 3. 

H 

Algorithm 3: Algorithm for finding a multiplicative inverse using simple operations. 

Input: a number 

Output: its multiplicative inverse 

INVS(z) 

(1) if z = throw an error. 

(2) Let x <— sign (z) , n <— 0. 

(3) while n < 50 

(4) Let x <- x (2 - zx) . 

(5) Let n <— n + 1. 



Example Problem 4.11. Devise a subroutine using only simple operations which computes the 
square root of an input number z. 

Solution: The temptation is to find a zero for f(x) = \fz — x. However, this equation is linear 
in x. Instead let f(x) = z — x 2 . You can easily see that if x is a positive root of f(x) = 0, then 
x = \fz. The Newton step becomes 



x k+ i <- x k 



4 



-2x k 

after some simplification this becomes 

x k , z 
Xk+1 ^ T 2x~ k 

Note this relies only on addition, multiplication and division. 
The final subroutine is given in Algorithm 4. 

H 



4.3 Secant Method 

The secant method for root finding is roughly based on Newton's method; however, it is not 
assumed that the derivative of the function is known, rather the derivative is approximated by 
the value of the function at some of the iterates, x k - More specifically, the slope of the tangent 
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Algorithm 4: Algorithm for finding a square root using simple operations. 
Input: a number 
Output: its square root 
SQRt(z) 

(1) if z < throw an error. 

(2) Let x <- 1, n <- 0. 

(3) while n < 50 

(4) Let x <- (a; + z/x) /2. 

(5) Let n n + 1. 




Figure 4.2: One iteration of the Secant method is shown for some quadratic function f(x). The 
secant line through (x k -i, f(xk~i)) and (xk,f(xk)) is shown. It happens to be the case that 
is smaller than \f(xk)\ , i.e., is a better guess than x^. 

line at (x^, f(xk)) , which is f'(xk) is approximated by the slope of the secant line passing through 
(x k -i,f(xk-i)) and (x k ,f(x k )) , which is 

f(Xk) ~ f(Xk-l) 
%k %k— 1 

Thus the iterate x/c+i is the root of this secant line. That is, it is a root to the equation 

— {x-x k ) = y- f{x k ). 
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Since the root has a y value of 0, we have 
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fjxk) ~ fi x k- 

Xk Xk—l 



(x k+1 - x k ) = -f(x k ), 

_ _ ( x k - Xk-1 
\f(xk) - f(x k - 



1, 



x k+ i = x k - 



Xk X k —i 



f(xk) - f(Xk-l) 



f(x k ), 
f{x k ). 



(4.2) 



You will note this is the recurrence of Newton's method, but with the slope f'(x k ) substituted 
by the slope of the secant line. Note also that the secant method requires two initial guesses, xo, x±, 
but can be used on a black box function. The secant method can suffer from some of the same 
problems that Newton's method does, as we will see. 

An iteration of the secant method is shown in Figure 4.2, along with the secant line. 



Example 4.12. Consider the secant method used on x 3 + x 2 — x — 1, with xo = 2, x\ 
Note that this function is continuous and has roots ±1. We give the iterates here: 



A: 


Xk 


f(xk) 





2 


9 


1 


0.5 


-1.125 


2 


0.666666666666667 


-0.925925925925926 


3 


1.44186046511628 


2.63467367653163 


4 


0.868254072087394 


-0.459842466254495 


5 


0.953491494113659 


-0.177482458876898 


6 


1.00706900811804 


0.0284762692197613 


7 


0.999661272951803 


-0.0013544492875992 


8 


0.999997617569723 


-9.52969840528617 x 10 - ' 


9 


1.0000000008072 


3.22880033820638 x 10~ 09 


10 


0.999999999999998 


-7.43849426498855 x 10 _ 



4.3.1 Problems 

As with Newton's method, convergence is dependent on f(x), and the initial estimates xq,X\. We 
illustrate a few possible problems: 



Example 4.13. Consider the secant method for the function f(x) = with x$ = 3, X\ = 4. 
As with Newton's method, the iterates diverge towards infinity, looking for a nonexistant root. We 
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rive some iterates here: 



k Xfc 

3 

1 4 

2 21.6548475770851 

3 33.9111765137635 

4 67.3380435135758 

5 117.820919458675 

6 210.543986613847 

7 366.889164762149 

8 637.060241341843 

9 1096.54125113444 

10 1878.34688714646 

11 3201.94672271613 

12 5437.69020766155 

13 9203.60222260594 

14 15533.1606791089 

15 26149.7196085218 

16 43924.8466075548 

17 73636.673898472 



f(xk) 

0.366204096222703 

0.346573590279973 

0.142011128224341 

0.103911011441661 

0.0625163004418104 

0.0404780904944712 

0.0254089165873003 

0.0160949419231219 

0.010135406045582 

0.00638363233543847 

0.00401318169994875 

0.0025208146648422 

0.00158175793894727 

0.000991714984152597 

0.000621298692241343 

0.000388975250950428 

0.000243375589137882 

0.000152191807070607 



Example 4.14. Consider Newton's method applied to the function f(x) = — — jj, with initial 
estimates xq = — l,xi = 1. 

You can easily verify that /(#o) = f(xi), and thus the secant line is horizontal. And thus X2 is not 
defined. 

Algorithm 5: Algorithm for finding root by secant method. 
Input: a function, initial guesses, an iteration limit, and a tolerance 
Output: a point for which the function has small value. 

run_secant(/, xO, xl, N, tol) 

(1) Let x <— xl, xp <— xO, fh <— /(xl), fp <— f(x0),n <— 0. 

(2) while n < N 

(3) if \fh\ < tol 

(4) return x. 

(5) Let fpx <- (fh - fp)/ (x - xp). 

(6) if \fpx\< tol 

(7) Warn "secant slope is too small; giving up." 

(8) return x. 

(9) Let xp <— x, fp <— //t. 

(10) Let x <— x — fh/ fpx. 

(11) Let fh<-f(x). ' 

(12) Letnf-n + 1. 



4.3.2 Convergence 

We consider convergence analysis as in Newton's method. We assume that r is a root of f(x), and 
let e ra = r — x n . Because the secant method involves two iterates, we assume that we will find some 
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relation between e k+ \ and the previous two errors, e k ,e k -i. 

Indeed, this is the case. Omitting all the nasty details (see Cheney & Kincaid [7]), we arrive at 
the imprecise equation: 

2/'(r) 

Again, the proof relies on finding some "attractor" region and using continuity. 

We now postulate that the error terms for the secant method follow some power law of the 
following type: 

|e fc+ i| ~ A\e k \ a . 

Recall that this held true for Newton's method, with a = 2. We try to find the a for the secant 
method. Note that 

II A I \C 

\ek\ = A\e k -i\ , 

so 

I I A-± I |i 

|efe_i| = A <* \e k \ a ■ 

Then we have 

A \e k \ a = |e fc+ i| = C \e k \ |e fc _i| = CA~* |e fe | 1+ ~ , 
Since this equation is to hold for all e k , we must have 

1 

a = 1 + -. 

a 

This is solved by a = \ (l + V5) w 1.62. Thus we say that the secant method enjoys superlinear 
convergence; This is somewhere between the convergence rates for bisection and Newton's method. 
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Exercises 

(4.1) Consider the bisection method applied to find the zero of the function f(x) = x 2 — 5x + 3, 
with ao = 0, &o = 1- What are a\, b\l What are 02, 62 ? 

(4.2) Approximate \/l0 by using two steps of Newton's method, with an initial estimate of xq = 3. 
(c/. Example Problem 4.11) Your answer should be correct to 5 decimal places. 

(4.3) Consider bisection for finding the root to cosx = 0. Let the initial interval Iq be [0, 2]. What 
is the next interval considered, call it I\! What is I-p. I%1 

(4.4) What does the sequence defined by 

1 1 

x = 1, Xfc+i = -x k H 

2 x k 

converge to? 

(4.5) Devise a subroutine using only simple operations that finds, via Newton's Method, the cubed 
root of some input number z. 

(4.6) Use Newton's Method to approximate \^9. Start with xq = 2. Find X2. 

(4.7) Use Newton's Method to devise a sequence xq,x\, . . . such that Xk — > In 10. Is this a reasonable 
way to write a subroutine that, given z, computes In z? (Hint: such a subroutine would require 
computation of e Xk . Is this possible for rational Xk without using a logarithm? Is it practical?) 

(4.8) Give an example (graphical or analytic) of a function, f(x) for which Newton's Method: 

(a) Does not find a root for some choices of xq. 

(b) Finds a root for every choice of xq. 

(c) Falls into a cycle for some choice of xq. 

(d) Converges slowly to the zero of f(x). 

(4.9) How will Newton's Method perform for finding the root to f(x) = y/\x\ = 0? 

(4.10) Implement the inverse finder described in Example Problem 4.10. Your m-file should have 
header line like: 

function c = invs(z) 

where z is the number to be inverted. You may wish to use the builtin function sign. As 
an extra termination condition, you should have the subroutine return the current iterate if 
zx is sufficiently close to 1, say within the interval (0.9999,0.0001). Can you find some z for 
which the algorithm performs poorly? 

(4.11) Implement the bisection method in octave/Matlab. Your m-file should have header line like: 
function c = run_bisection(f , a, b, tol) 

where f is the name of the function. Recall that feval(f ,x) when f is a string with the 
name of some function evaluates that function at x. This works for builtin functions and 
m-files. 

(a) Run your code on the function f(x) = cosx, with ao = 0, 60 = 2. In this case you can 
set f = "cos". 

(b) Run your code on the function f(x) = x 2 — 5x + 3, with ao = 0, 60 = 1- In this case you 
will have to set f to the name of an m-file (without the ".m") which will evaluate the 
given f(x). 

(c) Run your code on the function f(x) = x — cosx, with ao = 0, 60 = 1- 

(4.12) Implement Newton's Method. Your m-file should have header line like: 
function x = run_newton(f , fp, xO, N, tol) 

where f is the name of the function, and f p is its derivative. Run your code to find zeroes of 
the following functions: 
(a) f(x) = tanx — x. 
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(b) f(x) = x 2 — (2 + e)x + 1 + e, for e small. 

(c) f(x) = x 7 - 7x 6 + 21a; 5 - 35x 4 + 35x 3 - 21x 2 + 7x-l. 

(d) f(x) = (x-l) 7 . 

(4.13) Implement the secant method Your m-file should have header line like: 
function x = run_secant (f , xO, xl, N, tol) 

where f is the name of the function. Run your code to find zeroes of the following functions: 

(a) f(x)=t&nx — x. 

(b) f{x) = x 2 + 1. (Note: This function has no zeroes.) 

(c) f(x) = 2 + sinx. (Note: This function also has no zeroes.) 
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Chapter 5 

Interpolation 



5.1 Polynomial Interpolation 

We consider the problem of finding a polynomial that interpolates a given set of values: 



X 


x 






Xn 


y 


yo 


yi 




Vn 



where the Xi are all distinct. A polynomial p(x) is said to interpolate these data if p(xi) = yi for 
i = 0, 1, . . . , n. The x« values are called "nodes." 

Sometimes, we will consider a variant of this problem: we have some black box function, f(x), 
which we want to approximate with a polynomial p{x). We do this by finding the polynomial 
interpolant to the data 



X 


x 


Xl 




Xn 


m 


f(xo) 






f(Xn) 



for some choice of distinct nodes Xi. 



5.1.1 Lagranges Method 



As you might have guessed, for any such set of data, there is an n-degree polynomial that interpo- 
lates it. We present a constructive proof of this fact by use of Lagrange Polynomials. 

Definition 5.1 (Lagrange Polynomials). For a given set of n + 1 nodes Xj, the Lagrange 
polynomials are the n + 1 polynomials £i defined by 



(■i (Xj ) 



ifi^j 

1 if i = j 



Then we define the interpolating polynomial as 

n 

Pn{x) = ^2ydi(x). 



i=0 



If each Lagrange Polynomial is of degree at most n, then p n also has this property. The Lagrange 
Polynomials can be characterized as follows: 



ti(*)= n 



(5.1) 
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By evaluating this product for each Xj , we see that this is indeed a characterization of the Lagrange 
Polynomials. Moreover, each polynomial is clearly the product of n monomials, and thus has degree 
no greater than n. 

This gives the theorem 

Theorem 5.2 (Interpolant Existence and Uniqueness). Let {xj}™ =0 be distinct nodes. Then 
for any values at the nodes, {yi}™ =0 , there is exactly one polynomial, p(x) of degree no greater than 
n such that p(xi) = yi for i = 0, 1, . . . , n. 

Proof. The Lagrange Polynomial construction gives existence of such a polynomial p(x) of degree 
no greater than n. 

Suppose there were two such polynomials, call them p(x),q(x), each of degree no greater than 
n, both interpolating the data. Let r(x) = p(x) — q(x). Note that r(x) can have degree no greater 
than n, yet it has roots at xq, x\, . . . , x n . The only polynomial of degree < n that has n + 1 distinct 
roots is the zero polynomial, i.e., = r(x) = p(x) — q(x). Thus p(x),q(x) are equal everywhere, 
i.e., they are the same polynomial. □ 

Example Problem 5.3. Construct the polynomial interpolating the data 



X 


l 


l 

2 


3 


y 


3 


-10 


2 



by using Lagrange Polynomials. 
Solution: We construct the Lagrange Polynomials: 

e (x) 

h(x) 

Then the interpolating polynomial, in "Lagrange Form" is 

p 2 (x) = -3(x - l -)(x - 3) - 8(x - l)(x - 3) + 2 -{x - l)(x - \) 

H 

5.1.2 Newton's Method 

There is another way to prove existence. This method is also constructive, and leads to a different 
algorithm for constructing the interpolant. One way to view this construction is to imagine how 
one would update the Lagrangian form of an interpolant. That is, suppose some data were given, 
and the interpolating polynomial calculated using Lagrange Polynomials; then a new point was 
given (x n+ i, y n+ i), and an interpolant for the augmented set of data is to be found. Each Lagrange 
Polynomial would have to be updated. This could take a lot of calculation (especially if n is large). 

So the alternative method constructs the polynomials iteratively. Thus we create polynomials 
Pk(x) such that Pk{xi) = y% for < i < k. This is simple for k = 0, we simply let 



(x-i)(x-3) 



(1- 




-3) 


(x — 


l)(x 


-3) 


a- 




-3) 


(x — 


l)(x 





(3-l)(3-i) 



~(x ~ ^)(x - 3) 
|(* -!)(*-!) 



Po(x) = y , 
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-0.5 - 



0.5 - 



1.5 - 



-1 



- 



1 - 








0.5 



1.5 



2 



2.5 



3 



3.5 



Figure 5.1: The 3 Lagrange polynomials for the example are shown. Verify, graphically, that these 
polynomials have the property ii(xj) = 5ij for the nodes 1, |,3. 

the constant polynomial with value yo- Then assume we have a proper pk(x) and want to construct 
Pk+i(x). The following construction works: 



for some constant c. Note that the second term will be zero for any Xi for < i < k, so pk+i(x) 
will interpolate the data at xq, x±, . . . , x^. To get the value of the constant we calculate c such that 



Vk+i = Pk+i{xk+i) = Pk(xk+i) + c(xfc+i - x )(x k +i - xi) ■ ■ ■ - x k ). 



This construction is known as Newton's Algorithm, and the resultant form is Newton's form of 
the interpolant 

Example Problem 5.4. Construct the polynomial interpolating the data 



by using Newton's Algorithm. 

Solution: We construct the polynomial iteratively: 

p (x) = 3 

Pi(x) = 3 + c(x — 1) 
We want -10 = pi{\) = 3 + c(-|), and thus c = 26. Then 



Pk+i{x) = Pk(x) + c(x - x )(x - xx) ■ ■ ■ (x - x k ), 



x 1 A 3 



y 3 -10 2 
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0.5 1 1.5 2 2.5 3 3.5 4 



Figure 5.2: The interpolating polynomial for the example is shown. 

We want 2 = p 2 (3) = 3 + 26(2) + c(2)(§), and thus c = =j£. Then we get 

—53 1 
p 2 (x) = 3 + 26(x - 1) + -— (x - l)(x - -). 

5 2 

H 

Does Newton's Algorithm give a different polynomial? It is easy to show, by induction, that 
the degree of p n (x) is no greater than n. Then, by Theorem 5.2, it must be the same polynomial 
as the Lagrange interpolant. 

The two methods give the same interpolant, we may wonder "Which should we use?" Newton's 
Algorithm seems more flexible-it can deal with adding new data. There also happens to be a way 
of storing Newton's form of the interpolant that makes the polynomial simple to evaluate (in the 
sense of number of calculations required). 

5.1.3 Newton's Nested Form 

Recall the iterative construction of Newton's Form: 

Pk+i(x) = Pk(x) + c k (x - x )(x -xi)---(x- x k ). 
The previous iterate p k (x) was constructed similarly, so we can write: 

Vk+i{x) = \Pk-i{x) + c fc _i(x - x )(x - Xi) ■ ■ ■ (x - x k _i)] + c k (x - x )(x -xi)---(x- x k ). 
Continuing in this way we see that we can write 



n 






Pn(x) = y^ ck 






k=0 


0<j<k 





where an empty product has value 1 by convention. This can be rewritten in a funny form, where 
a monomial is factored out of each successive summand: 

p n (x) = C + (X - X )[CI + {X- Xl)[c 2 + (x - X 2 ) [. . .]]] 
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Supposing for an instant that the constants c& were known, this provides a better way of 
calculating p n (t) at arbitrary t. By "better" we mean requiring few multiplications and additions. 
This nested calculation is performed iteratively: 



V 1 = C n _i + {t- Xn^Vo 
V2 = Cn-2 + (t - X n -2)Vl 

V n = Co + (t- x )v n -i 

This requires only n multiplications and 2n additions. Compare this with the number required 
for using the Lagrange form: at least n 2 additions and multiplications. 



5.1.4 Divided Differences 

It turns out that the coefficients Ck for Newton's nested form can be calculated relatively easily by 
using divided differences. We assume, for the remainder of this section, that we are considering 
interpolating a function, that is, we have values of f(xi) at the nodes Xj. 

Definition 5.5 (Divided Differences). For a given collection of nodes {xj}" =0 and values 
{f( x i)}x=oi a ^ th order divided difference is a function of k + 1 (not necessarily distinct) nodes, 
written as 

/ [ x ioi x ii > • • • > x ik\ 

The divided differences are defined recursively as follows: 

• The th order divided differences are simply defined: 

f[Xi] = f(Xi). 

• Higher order divided differences are the ratio of differences: 

p r I / [ x ii i x i2 i • • • i x i}~\ ~ f \_ x io i x ii ) • • • ) 

/ [ x io ) x h > ■ ■ ■ > x ik\ 



We care about divided differences because coefficients for the Newton nested form are divided 
differences: 



Ck = f [x ,xi,.. .,x k \. 



(5.3) 



Because we are only interested in the Newton method coefficients, we will only consider divided 
differences with successive nodes, i.e., those of the form / [xj, Xj+i, . . . , Xj+fc]. In this case the 
higher order differences can more simply be written as 



/ [Xj, Xj+l, . . . , Xj+k] — 



f [ x j+l, x j+2, x j+k] ~ f [Xj,Xj + l, ... , X j+k _{\ 

x j+k x j 
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The graphical way to calculate these things is a "pyramid scheme" 1 , where you compute the 
following table columnwise, from left to right: 



X 


f[] 


/[,] 


/[,,] 


/[,,,] 


x 


f[xo] 














f[xo,Xi] 








X\ 


f[xi] 




f[xo,xi 


x 2 ] 








f[xi,x 2 ] 






f[x ,Xl,X 2 ,X3,] 


X2 


f[x2] 




f[xi,x 2 


X3] 








f[x2,X 2 ] 








X3 


/N 











Note that by definition, the first column just echoes the data: f[xi] = f{xi). 
We drag out our example one last time: 

Example Problem 5.6. Find the divided differences for the following data 



X 


1 


1 

2 


3 


Si?) 


3 


-10 


2 



Solution: We start by writing in the data: 



X 


/[] 


/[,] 


/[,,] 


1 


3 






1 

2 


-10 






3 


2 







Then we calculate: 

f[xo,xi 



(1/2) - 1 
Adding these to the table, we have 



2d, and j[xi,x 2 \ 



3 - (1/2) 



24/5. 



X 


/[] 


/[,] 


1 


3 








26 


1 

2 


-10 








24 






5 


3 


2 





/[,, 



Then we calculate 



So we complete the table: 



f[x ,xi,x 2 ] 



24/5 - 26 _ -53 
3-1 ~ ~~5~ ' 



X 


/[] 


/[,] 


/[,,] 


1 


3 










26 




1 

2 


-10 




-53 

5 






24 








5 




3 


2 







1 That's supposed to be a joke. 
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You should verify that along the top line of this pyramid you can read off the coefficients for 
Newton's form, as found in Example Problem 5.4. H 

5.2 Errors in Polynomial Interpolation 

We now consider two questions: 

1. If we want to interpolate some function f(x) at n + 1 nodes over some closed interval, how 
should we pick the nodes? 

2. How accurate can we make a polynomial interpolant over a closed interval? 

You may be surprised to find that the answer to the first question is not that we should make 
the Xi equally spaced over the closed interval, as the following example illustrates. 

Example 5.7 (Runge Function). Let 

f(x) = {l + x 2 )-\ 

(known as the Runge Function), and let p n (x) interpolate / on n equally spaced nodes, including 
the endpoints, on [—5,5]. Then 

lim max \p n (x) — f(x)\ = oo. 

n->oo x e[-5,5] 

This behaviour is shown in Figure 5.3. 

It turns out that a much better choice is related to the Chebyshev Polynomials ("of the first 
kind"). If our closed interval is [—1,1], then we want to define our nodes as 



Xi = cos 



2i + l 
2n + 2 



7T 



< i < n. 



Literally interpreted, these Chebyshev Nodes are the projections of points uniformly spaced on 
a semi circle; see Figure 5.4. By using the Chebyshev nodes, a good polynomial interpolant of the 
Runge function of Example 5.7 can be found, as shown in Figure 5.5. 

5.2.1 Interpolation Error Theorem 

We did not just invent the Chebyshev nodes. The fact that they are "good" for interpolation 
follows from the following theorem: 

Theorem 5.8 (Interpolation Error Theorem). Let p be the polynomial of degree at most n 
interpolating function / at the n+1 distinct nodes xq, xi, . . . , x n on [a, b]. Let /( n+1 ) be continuous. 
Then for each x £ [a, b] there is some £ G [a, b] such that 



1 " 
/(*) " V{x) = 1 ——fW (0 H (x - Xl ) . 



(n + 1)! 

You should be thinking that the term on the right hand side resembles the error term in Taylor's 
Theorem. 
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0.8 



0.6 



0.4 



0.2 



-0.2 



-0.4 



runge 
3 nodes 
7 nodes 

10 nodes 




-2 2 

(a) 3,7,10 equally spaced nodes 




(b) 15 equally spaced nodes 



(c) 50 equally spaced nodes 



Figure 5.3: The Runge function is poorly interpolated at n equally spaced nodes, as n gets large. 
At first, the interpolations appear to improve with larger n, as in (a). In (b) it becomes apparent 
that the interpolant is good in center of the interval, but bad at the edges. As shown in (c), this 
gets worse as n gets larger. 
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Figure 5.4: The Chebyshev nodes are the projections of nodes equally spaced on a semi circle. 




(a) 3,7,10 Chebyshev nodes (b) 17 Chebyshev nodes 



Figure 5.5: The Chebyshev nodes yield good polynomial interpolants of the Runge function. Com- 
pare these figures to those of Figure 5.3. For more than about 25 nodes, the interpolant is indis- 
tinguishable from the Runge function by eye. 

Proof. First consider when x is one of the nodes xf, in this case both the LHS and RHS are zero. 
So assume x is not a node. Make the following definitions 

n 

w{t) = n (t - , 

= f{x) -p(x) 
w(x) 

0(t) = f(t)-p(t)-cw(t). 

Since x is not a node, w(x) is nonzero. Now note that 4>(xi) is zero for each node Xi, and that by 
definition of c, that 4>(x) = for our x. That is 4>(t) has n + 2 roots. 

Some other facts: f,p,w have n + 1 continuous derivatives, by assumption and definition; thus 
cf) has this many continuous derivatives. Apply Rolle's Theorem to cj)(t) to find that (j)'(t) has n + 1 
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roots. Then apply Rolle's Theorem again to find 4>"(t) has n roots. In this way we see that (f) (n+l \t) 
has a root, call it £. That is 

o = 4> {n+1) (0 = f (n+1) (0-p {n+1) (0 - 

But p(t) is a polynomial of degree < n, so j/ ra+1 ) is identically zero. And w(t) is a polynomial 
of degree n + 1 in t, so its n + 1 th derivative is easily seen to be (n + 1)! Thus 

= /(" +1 )(e)-c(n + l)! 
c(n + l)! = /(" +1 )(e) 



/(a:) -p(x) 
w(x) 



1 



/(»+!)(£) 



(n + 1)!' 

which is what was to be proven. 

Thus the error in a polynomial interpolation is given as 



□ 



1 n 
/(*) " V{x) = j—-^fW (£) H (x - x t ) . 



We have no control over the function f(x) or its derivatives, and once the nodes and / are fixed, p 
is determined; thus the only way we can make the error \f(x) — p(x)\ small is by judicious choice 
of the nodes Xj. 

The Chebyshev nodes on [—1, 1] have the remarkable property that 



^=0 



< 2" 



for any t £ [—1, 1] . Moreover, it can be shown that for any choice of nodes Xi that 



max 
*e[-i,i] 



II<« 



> 2~ n . 



Thus the Chebyshev nodes are considered the best for polynomial interpolation. 

Merging this result with Theorem 5.8, the error for polynomial interpolants defined on Cheby- 
shev nodes can be bounded as 

I/O) ~p(x)\ < 



max 



2 n (n + l)Ue[-i,i] 



f (n+1 Ho 



The Chebyshev nodes can be rescaled and shifted for use on the general interval [a, (3\ . In this 
case they take the form 



(3 — a 
Xi = cos 



2i + l 
2n + 2 



7T 



a + (3 



0<i<n. 



In this case, the rescaling of the nodes changes the bound on T\(t — %i) so the overall error bound 
becomes 

08 - «) n+1 



\f(x)-p(x)\ < 



max 



2 2n+l ( n+ l)! 5e[Qj/3] 



f(n + l) {0 



for x £ [a, (3] . 
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Example Problem 5.9. How many Chebyshev nodes are required to interpolate the function 
f(x) = sin(x) + cos(x) to within 10~ 8 on the interval [0, 7r]? 
Solution: We first find the derivatives of / 

f'(x) = cos(x) — sin(x) 
/ (x) = — sin(x) — cos(x) 
f"'(x) = — cos(x) + sin(x) 



As a crude approximation we can assert that 



< Icosxl + I sin x I < 2. 



Thus it suffices to take n large enough such that 



:2 < 10- 8 



2 2n+l ( n+ 1)!' 

By trial and error we see that n = 10 suffices. H 
5.2.2 Interpolation Error for Equally Spaced Nodes 

Despite the proven superiority of Chebyshev Nodes, and the problems with the Runge Function, 
equally spaced nodes are frequently used for interpolation, since they are easy to calculate 2 . We 
now consider bounding 

n 

max TT |x — Xjl , 

where 

Xi = a + hi = a + — — —i, i = 0, 1, . . . , n. 

n 

Start by picking an x. We can assume x is not one of the nodes, otherwise the product in 
question is zero. Then x is between some Xj, Xj + \ We can show that 

< — . 

by simple calculus. 

Now we claim that \x — Xj| < (j — i + 1) h for i < j, and |x — x«| < (i — j) h for j + 1 < i. Then 

n h 2 

Y[\x-Xi\< — [(j + 1)W] [(n - jW'i- 1 ] . 

It can be shown that (j + l)!(n — j)\ < n\, and so we get an overall bound 

h n+1 n\ 

X — X{\ < 

i=0 



Y\_ \X - Xi\ < 



2 Never underestimate the power of laziness. 



74 



CHAPTER 5. INTERPOLATION 



The interpolation theorem then gives us 



h n+1 

\f(x) —p(x)\ < —, — max 



f(n + l) {0 



where h = (b — a)/n. 

The reason this result does not seem to apply to Runge's Function is that for Runge's 
Function becomes unbounded as n — > oo. 

Example Problem 5.10. How many equally spaced nodes are required to interpolate the function 

f(x) = sin(x) + cos(a;) to within 10 -8 on the interval [0, 7r]? 

Solution: As in Example Problem 5.9, we make the crude approximation 



/<*>(*) 



< Icosxl + I sin a; I < 2. 



Thus we want to make n sufficiently large such that 

h n+i 



-2 < 10"°, 



4(n+ 1) 

where h = (tt — 0)/n. That is we want to find n large enough such that 

7T ri+1 



2n n+1 (n + 1) 

By simply trying small numbers, we can see this is satisfied if n = 12. 
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Exercises 

(5.1) Find the Lagrange Polynomials for the nodes { — 1, 1}. 

(5.2) Find the Lagrange Polynomials for the nodes { — 1, 1, 5}. 

(5.3) Find the polynomial of degree no greater than 2 that interpolates 



X 


-l 


1 


5 


y 


3 


3 


-2 



(5.4) Complete the divided differences table: 



X 


/[] 


/[,] 


/[,,] 


-1 


3 






1 


3 






5 


-2 







Find the Newton form of the polynomial interpolant. 
(5.5) Find the polynomial of degree no greater than 3 that interpolates 



X 


-l 


1 


5 


-3 


y 


3 


3 


-2 


4 



(Hint: reuse the Newton form of the polynomial from the previous question.) 
(5.6) Find the nested form of the polynomial interpolant of the data 



X 


l 


3 


4 


6 


y 


-3 


13 


21 


1 



by completing the following divided differences table: 



X 


/[] 


/[,] 


/[,,] 


/[,,,] 


1 


-3 








3 


13 








4 


21 








6 


1 









(5.7) Find the polynomial of degree no greater than 3 that interpolates 



X 


l 





3/2 


2 


y 


3 


2 


37/8 


8 



(5.8) Complete the divided differences table: 



X 


/[] 


/[,] 


/[,,] 


/[,,,] 


-1 


6 











3 








1 


2 








2 


3 
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(Something a little odd should have happened in the last column.) Find the Newton form of 
the polynomial interpolant. Of what degree is the polynomial interpolant? 

(5.9) Let p(x) interpolate the function cos(x) at n equally spaced nodes on the interval [0, 2] . Bound 
the error 

max \p(x) — cos(x)| 

0<x<2 

as a function of n. How small is the error when n = 10? How small would the error be if 
Chebyshev nodes were used instead? How about when n = 10? 

(5.10) Let p(x) interpolate the function x~ 2 at n equally spaced nodes on the interval [0.5,1]. 
Bound the error 

max \p(x) — x~ 2 \ 

0.5<a;<l 1 V ' 1 

as a function of n. How small is the error when n = 10? 

(5.11) How many Chebyshev nodes are required to interpolate the function ^ to within 10~ 6 on 
the interval [1,2]? 

(5.12) Write code to calculate the Newton form coefficients, by divided differences, for the nodes 
Xi and values /(xj). Your m-file should have header line like: 

function coefs = newtonCoef (xs, fxs) 

where xs is the vector of n + 1 nodes, and f xs the vector of n + 1 values. Test your code on 
the following input: 

octave :1> xs = [1-12-23-34 -4] ; 
octave :2> fxs = [1 1 2 3 5 8 13 21] ; 
octave :3> newtonCoef (xs, fxs) 
ans = 



1.00000 -0.00000 0.33333 -0.08333 0.05000 0.00417 -0.00020 -0.00040 



(a) What do you get when you try the following? 
octave :4> xs = [1 4 5 9 14 23 37 60]; 
octave :5> fxs = [3 141592 6]; 
octave :6> newtonCoef (xs ,ys) 

(b) Try the following: 

octave :7> xs = [1 3 4 2 8 -2 14 23 15]; 
octave :8> fxs = xs.*xs + xs .+4; 
octave :9> newtonCoef (xs, fxs) 
(5.13) Write code to calculate a polynomial interpolant from its Newton form coefficients and the 
node values. Your m-file should have header line like: 
function y = calcNewton(t , coefs ,xs) 

where coefs is a vector of the Newton coefficients, xs is a vector of the nodes Xj, and y is 
the value of the interpolating polynomial at t. Check your code against the following values: 
octave :1> xs = [13 4]; 
octave :2> coefs = [6 5 1]; 
octave :3> calcNewton (5,coefs,xs) 
ans = 34 

octave :4> calcNewton (-3 , coef s ,xs) 
ans = 10 

(a) What do you get when you try the following? 
octave :5> xs= [3 14592687 0]; 
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octave :6> coefs = [1-12-23-34-45 -5]; 
octave : 7> calcNewton (0 . 5 , coefs , xs) 
(b) Try the following 

octave :8> xs= [3 14592687 0]; 

octave :9> coefs = [1-12-23-34-45 -5]; 

octave : 10> calcNewton(l , coef s ,xs) 
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Chapter 6 

Spline Interpolation 



Splines are used to approximate complex functions and shapes. A spline is a function consisting of 
simple functions glued together. In this way a spline is different from a polynomial interpolation, 
which consists of a single well defined function that approximates a given shape; splines are normally 
piecewise polynomial. 



6.1 First and Second Degree Splines 



Splines make use of partitions, which are a way of cutting an interval into a number of subintervals. 

Definition 6.1 (Partition). A partition of the interval [a, b] is an ordered sequence {tj}™ =0 such 
that 

a = t < h < ■ ■ ■ < t n -i <t n = b 

The numbers ti are known as knots. 

A spline of degree 1, also known as a linear spline, is a function which is linear on each subinterval 
defined by a partition: 

Definition 6.2 (Linear Splines). A function S is a spline of degree 1 on [a, b] if 

1. The domain of S is [a, b\. 

2. S is continuous on [a, b\. 

3. There is a partition {ii}™ =0 of [a,b] such that on each [ti, , S is a linear polynomial. 
A linear spline is defined entirely by its value at the knots. That is, given 



t 


to 


h 






y 


yo 


yi 




y n 



there is only one linear spline with these values at the knots and linear on each given subinterval. 
For a spline with this data, the linear polynomial on each subinterval is defined as 

Vi+i - yi 



Si(x) = yi + 



U+l ~ ti 



(x - ti) 



Note that if x € [U, U+i] , then x — ti > 0, but x — tj_i < 0. Thus if we wish to evaluate S(x), we 
search for the largest i such that x — ti > 0, then evaluate Si(x). 

Example 6.3. The linear spline for the following data 



t 


0.0 


0.1 


0.4 


0.5 


0.75 


1.0 


y 


1.3 


4.5 


2.0 


2.1 


5.0 


3 



is shown in Figure 6.1. 
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6 



5 - 



4 - 



3 - 



2 - 




1 1 1 1 1 1 1 1 

0.2 0.4 0.6 0.8 1 

Figure 6.1: A linear spline. The spline is piecewise linear, and is linear between each knot. 
6.1.1 First Degree Spline Accuracy 

As with polynomial functions, splines are used to interpolate tabulated data as well as functions. In 
the latter case, if the spline is being used to interpolate the function /, say, then this is equivalent 
to interpolating the data 



t 


to 


h 




tn 


V 


f(to) 






f(tn) 



A function and its linear spline interpolant are shown in Figure 6.2. The spline interpolant in 
that figure is fairly close to the function over some of the interval in question, but it also deviates 
greatly from the function at other points of the interval. We are interested in finding bounds on 
the possible error between a function and its spline interpolant. 

To find the error bound, we will consider the error on a single interval of the partition, and use 
a little calculus. 1 Suppose p(t) is the linear polynomial interpolating f(t) at the endpoints of the 
subinterval then for t G [ij,tj + i] , 

\f(t) - p(t)\ < max{|/(i) - f(U)\ , \f(t) - f(t i+ i)\} . 

That is, \f(t) — p(t)\ is no larger than the "maximum variation" of f(t) on this interval. 
In particular, if f'(t) exists and is bounded by M\ on [if, tj+i], then 

\f(t)-p(t)\<^(t i+1 -U). 

Similarly, if f"(t) exists and is bounded by M2 on then 

\f(t)-p(t)\<^(t l+1 - ti ) 2 . 

Over a given partition, these become, respectively 



Although we could directly claim Theorem 5.8, it is a bit of overkill. 
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Figure 6.2: The function f(t) = 4.1 + sin (1/(0. 08i + 0.5)) is shown, along with the linear spline 
interpolant, for the knots 0,0.1,0.4,0.6,0.75,0.9,1.0 For some values of t, the function and its 
interpolant are very close, while they vary greatly in the middle of the interval. 





-p(t)\ < 


Mi 

~Y 


max (ti+i -U) , 

0<i<n 


\f(t) 


-v{t)\ < 


M 2 
8 


max (t i+1 - ti) 2 . 

0<i<n 







(6.1) 



If equally spaced nodes are used, these bounds guarantee that spline interpolants become better 
as the number of nodes is increased. This contrasts with polynomial interpolants, which may get 
worse as the number of nodes is increased, cf. Example 5.7. 

6.1.2 Second Degree Splines 

Piecewise quadratic splines, or splines of degree 2, are defined similarly: 

Definition 6.4 (Quadratic Splines). A function Q is a quadratic spline on [a,b] if 

1. The domain of Q is [a,b]. 

2. Q is continuous on [a, b]. 

3. Q' is continuous on (a,b). 

4. There is a partition {ti}" =0 of [a, b] such that on Q is a polynomial of degree at most 
2. 



Example 6.5. The following is a quadratic spine: 

Q(x) = 



x 



—x x < 0, 
2 x < x < 2, 



-x 2 + 7x-8 2<x. 

Unlike linear splines, quadratic splines are not defined entirely by their values at the knots. We 
consider why that is. The spline Q{x) is defined by its piecewise polynomials, 

Qi(x) = diX 2 + biX + Cj. 
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Thus there are 3n parameters to define Q(x). 
For each of the n subintervals, the data 



t 


to 


h 






V 


yo 


yi 




Vn 



give two equations regarding Qi(x), namely that Qi(U) must equal yi and must equal j/i+i- 

This is 2n equations. The condition on continuity of Q' gives a single equation for each of the n — 1 
internal nodes. This totals 3n — 1 equations, but 3n unknowns. This system is under determined. 

Thus some additional user-chosen condition is required to determine the quadratic spline. One 
might choose, for example, Q'(a) = 0, or Q"(a) = 0, or some other condition. 

6.1.3 Computing Second Degree Splines 

Suppose the data 



t 


to 


h 




tn 


y 


yo 


yi 




y n 



are given. Let Zi = Q'^U), and suppose that the additional condition to define the quadratic spline 
is given by specifying zq. We want to be able to compute the form of Qi(x). 
Because Qi(U) = yi, Q'^U) = Zi, Q'^U+i) = we see that we can define 

Qi( x ) = wrr- — tt ( x ~ + z i(x~ U) + yi. 

* {H+1 — H) 

Use this at tj+i : 

y i+ i = Qi(t i+ i) = — ^— (t i+ i - ti) 2 + Zi (t i+ i - U) + y i} 

* {H+l — ti) 

yi+i -yi = Zt+1 2 (U+i - ti) + Zi (t i+ i - u) , 

Vi+i ~Vi = 2 ~ tl > ■ 

Thus we can determine, from the data alone, Zi + \ from Zf. 

_ n yi+i - yi 

Zi+l — Zi. 

H+l — ti 

6.2 (Natural) Cubic Splines 

If you recall the definition of the linear and quadratic splines, probably you can guess the definition 
of the spline of degree k: 

Definition 6.6 (Splines of Degree k). A function S is a spline of degree k on [a, b] if 

1. The domain of S is [a, b\. 

2. S, S', S", S^- 1 ) are continuous on (a, b). 

3. There is a partition {tj}™ =0 of [a, b] such that on [ti, U+i], S is a polynomial of degree < k. 

You would also expect that a spline of degree k has k — 1 "degrees of freedom," as we show 
here. If the partition has n + 1 knots, the spline of degree k is defined by n(k + 1) parameters. The 
given data 



t 


to 


h 




tn 


y 


yo 


yi 




y n 
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provide 2n equations. The continuity of S', S" , . . . , S^^ 1 ^ at the n — 1 internal knots gives (k — 
l)(n — 1) equations. This is a total of n(k + 1) — (k — 1) equations. Thus we have k — 1 more 
unknowns than equations. Thus, barring some singularity, we can (and must) add k — 1 constraints 
to uniquely define the spline. These are the degrees of freedom. 

Often k is chosen as 3. This yields cubic splines. We must add 2 extra constraints to define the 
spline. The usual choice is to make 

S"(to) = S"(t n ) = 0. 

This yields the natural cubic spline. 



6.2.1 Why Natural Cubic Splines? 

It turns out that natural cubic splines are a good choice in the sense that they are the "inter- 
polant of minimal H 2 seminorm." The corollary following this theorem states this in more easily 
understandable terms: 

Theorem 6.7. Suppose / has two continuous derivatives, and S is the natural cubic spline inter- 
polating / at knots a = to < t\ < . . . < t n = b. Then 



f [S"{x)] 2 dx < f [f"{x)] 2 dx 

J a J a 



Proof. We let g{x) = f(x) — S(x). Then g(x) is zero on the (n + 1) knots t{. Derivatives are linear, 
meaning that 

f"(x) = S"(x)+g"(x). 

Then b b b b 

f [f"(x)] 2 dx= f [S"(x)] 2 dx+ f [g"(x)] 2 dx+ f 2S"(x)g"(x)dx. 

J a J a J a J a 

We show that the last integral is zero. Integrating by parts we get 
f S"(x)g"(x)dx = S"g' 

J a 

because S"(a) = S"(b) = 0. Then notice that S is a polynomial of degree < 3 on each interval, thus 
S"'(x) is a piecewise constant function, taking value Cj on each interval [ij,tj + i]. Thus 

n— 1 „h n—1 



r-b b pb rb 

S"'g'dx = - / S'"g'dx, 



/ s'" g ' dx = Y, c l9 ' dx = ag\ = o, 



i=0 " u i=0 

with the last equality holding because g(x) is zero at the knots. □ 

Corollary 6.8. The natural cubic spline is best twice-continuously differentiable interpolant for a 
twice-continuously differentiable function, under the measure given by the theorem. 

Proof. Let / be twice-continuously differentiable, and let S be the natural cubic spline interpolating 
f(x) at some given nodes {ij}" =0 . Let R(x) be some twice-continuously differentiable function 
which also interpolates f(x) at these nodes. Then S(x) interpolates R(x) at these nodes. Apply 
the theorem to get 



f [S"(x)} 2 dx< f [R"{x)] 2 dx 

J a J a 



□ 
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6.2.2 Computing Cubic Splines 

First we present an example of computing the natural cubic spline by hand: 
Example Problem 6.9. Construct the natural cubic spline for the following data: 



t 


-l 





2 


y 


3 


-1 


3 



Solution: The natural cubic spline is defined by eight parameters: 

S(x) = 

We interpolate to find that d = h = — 1 and 



ax 3 + bx 2 + cx + d x G [-1, 0] 
ex 3 + fx 2 + gx + h x G [0, 2] 



We take the derivative of S: 



S'(x) = 



—a + 6 — c — 1 = 3 
8e + 4/ + 2g - 1 = 3 



3ax 2 + 2bx + c x G [-1,0] 
3ex 2 + 2fx + g x G [0, 2] 



Continuity at the middle node gives c = g. Now take the second derivative of 5: 

„ f 6ax + 2b x G [-1,0] 

lX j ~ \ 6ex + 2/ x G [0, 2] 

Continuity at the middle node gives b = f. The natural cubic spline condition gives —6a + 2b = 
and 12e + 2f = 0. Solving this by "divide and conquer" gives 

o/ n J ar 3 + 3x 2 - 2x - 1 xG[-l,0] 
(Xj ~ \ -\x 3 + 3x 2 - 2x - 1 x G [0, 2] 

H 

Finding the constants for the previous example was fairly tedious. And this is for the case of 
only three nodes. We would like a method easier than setting up the \n equations and unknowns, 
something akin to the description in Subsection 6.1.3. The method is rather tedious, so we leave it 
to the exercises. 



6.3 B Splines 

The B splines form a basis for spline functions, whence the name. We presuppose the existence of 
an infinite number of knots: 

. . . < t 2 < h < t < h < t 2 < ■ ■ ■ , 

with lim^-oo tk = —oo and lim^oo tk = oo. 

The B splines of degree are defined as single "blocks" : 



S?(x) 



1 ti < X < t i+ l 

otherwise 
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The zero degree B splines are continuous from the right, are nonzero only on one subinterval 
[U,ti + i), sum to 1 everywhere. 

We justify the description of B splines as basis splines: If S is a spline of degree on the given 
knots and is continuous from the right then 



That is, the basis splines work in the same way that Lagrange Polynomials worked for polynomial 
interpolation. 

The B splines of degree k are defined recursively: 



Some B splines are shown in Figure 6.3. 

The B splines quickly become unwieldy. We focus on the case k = 1. The B spline Bj(x) is 

• Piecewise linear. 

• Continuous. 

• Nonzero only on (^,£3+2). 



These B splines are sometimes called hat functions. Imagine wearing a hat shaped like this! What- 
ever. 

The nice thing about the hat functions is they allow us to use analogy. Harken back to polyno- 
mial interpolation and the Lagrange Functions. The hat functions play a similar role because 



S(x) = Y,S(x l )B?(x). 




• 1 at ij+i. 




Then if we want to interpolate the following data with splines of degree 1: 



t to t\ 



t 



n 



v yo vi 



We can immediately set 



n 



s(x) = Y,y l Bl 1 {x). 



i=0 
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Figure 6.3: Some B splines for the knots to = 1, ti = 2.5, £2 = 3.5,^3 = 4 are shown. In (a), one of 
the degree B splines; in (b), a degree 1 B spline; in (c), two of the degree 1 B splines and a degree 
2 B spline are shown. The two hat functions "merge" together to give the quadratic B spline. 
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Exercises 



S(x) = 



S(x) 



on [0,4]? 


Why or why 


~3x + 2 


: < x < 1 


-2x + 4 


: 1 < x < 4 


on [0,2]? 


Why or why 


fx + 3 


: < x < 1 


I 3 


: 1 < x < 2 


on [0,4]? 


Why or why 


fx 2 + 3 


: < x < 3 


1 5x — 6 


: 3 < x < 4 



(6.4) Find constants, a, [3 such that the following is a linear spline on [0,5]. 



S{x) 



Ax -2 
ax + (5 
-2x + 10 



< x < 1 

1 < x < 3 
3 < x < 5 



(6.5) Is the following function a quadratic spline on [0, 4]? Why or why not? 

, fx 2 + 3 : < x < 3 
Q(x) = < 

[5x-6 : 3 < x < 4 

(6.6) Is the following function a quadratic spline on [0, 2]? Why or why not? 

'x 2 + 3x + 2 :0<x<l 



Q(x) 



2x 2 + x + 3 :l<x<2 



(6.7) Find constants, a, (3, 7 such that the following is a quadratic spline on [0,5]. 

Q(x) = < 

(6.8) Find the quadratic spline that interpolates the following data: 



'ix 2 + 2x + § 
ax 2 + 0x + 7 
3x 2 - 7x + 12 



< x < 1 

1 < x < 3 
3 < x < 5 



t 





1 


4 


y 


1 


-2 


1 



To resolve the single degree of freedom, assume that Q'(0) = —Q'(4). Assume your solution 
takes the form 

|ai(x-l) 2 + /3i(x-l)-2 :0<x<l 
|a 2 (x-l) 2 + /3 2 (x-l)-2 :l<x<4 

Find the constants a±, Pi, a 2 , /?2- 



Q(x) 
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(6.9) Find the natural cubic spline that interpolates the data 



X 





1 


3 


y 


4 


2 


7 



It may help to assume your answer has the form 

S(X): 



Ax 3 + Bx 2 + Cx + 4 
D(x - l) 3 + E(x - l) 2 + F(x-l) + 2 



< x < 1 
1< x < 3 



Chapter 7 

Approximating Derivatives 



7.1 Finite Differences 

Suppose we have some blackbox function f(x) and we wish to calculate f'{x) at some given x. Not 
surprisingly, we start with Taylor's theorem: 

f( x + h ) = f(x) + f'(x)h+^-. 

Rearranging we get 

_ f( x + h)-f(x) _ f"(Qh 
nX) ~ h 2 ' 

Remember that £ is between x and x + h, but its exact value is not known. If we wish to calculate 
f'(x), we cannot evaluate /"(£), so we approximate the derivative by dropping the last term. That 
is, we calculate [f(x + h) — f(x)] /h as an approximation 1 to f'(x). In so doing, we have dropped 
the last term. If there is a finite bound on f"(z) on the interval in question then the dropped term 
is bounded by a constant times h. That is, 



fix) = + + O (h) 



(7.1) 



The error that we incur when we approximate f'(x) by calculating [f(x + h) — f(x)] /his called 
truncation error. It has nothing to do with the kind of error that you get when you do calculations 
with a computer with limited precision; even if you worked in infinite precision, you would still 
have truncation error. 

The truncation error can be made small by making h small. However, as h gets smaller, precision 
will be lost in equation 7.1 due to subtractive cancellation. The error in calculation for small h 
is called roundoff error. Generally the roundoff error will increase as h decreases. Thus there is 
a nonzero h for which the sum of these two errors is minimized. See Example Problem 7.5 for an 
example of this. 

The truncation error for this approximation is O (h). We may want a more precise approxima- 
tion. By now, you should know that any calculation starts with Taylor's Theorem: 

f(x + h) = f( x ) + f'( x )h + ^lh 2 + ^^h 3 

f(x-h) = f(x)-f'(x)h + ^lh 2 -^Mh" 



lr This approximation for f'(x) should remind you of the definition of f'(x) as a limit. 
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By subtracting these two lines, we get 

f{x + h)- f(x - h) 

Thus 

2f'(x)h 
f'(x) 



2f\x)h+ f " { ^ )+ J" {i2) h Z . 



f(x + h)- f{x - h) 

fjx + h)- fjx - h) 
2h 



r«i)+r(6) fc s 



3! 



h 2 



f ,n (P \4-f'"(P \ 

If f (x) is continuous, then there is some f between £ 1; £ 2 such that /"'(£) = 1 ^ 2> . (This 

is the MVT at work.) Assuming some uniform bound on /'"(■), we get 



/ , W = /( ' + lk) - /(x -* ) +o(tf) 



(7.2) 



In some situations it may be necessary to use evaluations of a function at "odd" places to 
approximate a derivative. These are usually straightforward to derive, involving the use of Taylor's 
Theorem. The following examples illustrate: 

Example Problem 7.1. Use evaluations of / at x + h and x + 2h to approximate fix), assuming 
fix) is an analytic function, i.e., one with infintely many derivatives. 

Solution: First use Taylor's Theorem to expand f{x + h) and f{x + 2h), then subtract to get some 
factor of fix): 



fix + 2h) 
fix + h) 



fix) + 2hfix) + yg-fix) + ^f'ix) + ^f 4 \x 



+ 



fix) + hfix) + % f\x) + %f"ix) + %f 4 \x) + ... 



fix + 2h)-fix + h) = hfix) + ^fix)+ 7 -^f'ix) + 1 -^f 4 \x) + . 
ifix + 2h)-fix + h))/h = f(x) + fr(x) + ^r(x) + ifi/( 4 )(x) + ... 

Thus ifix + 2h) - fix + h)) /h = fix) + O ih) 
Example Problem 7.2. Show that 

4fix + h)-fix + 2h)-3fix) 



2h 



= fix) + 0(h 2 ), 



for / with sufficient number of derivatives 

Solution: In this case we do not have to find the approximation scheme, it is given to us. We only 
have to expand the appropriate terms with Taylor's Theorem. As before: 

fix + h) = /(x) + TO + f/"(x) + S/'"(x) + ... 
Afix + h) = Afix) + Ahfix) + ^f"ix) + ^f"ix) + ... 
fix + 2h) = fix) + 2hfix) + ^fix) + ^nx) + ... 



4/(s + h)- fix + 2h) 
Afix + h)- fix + 2h) - 3/(x) 
i4fix + h)-fix + 2h)-3fix))/2h 



3/(x) + 2hfix) + Ofix) + =^f"ix) + 



2hf'[x) + + . . 



f'(x) + =^-f'"(x) + ... 
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7.1.1 Approximating the Second Derivative 

Suppose we want to approximate the second derivative of some blackbox function f(x). Again, 
start with Taylor's Theorem: 

f(x + h) = f(x) + f\x)h+^h 2 + ^h 3 + ^h 4 + ... 
f(x-h) = f(x)-f\x)h + l^hi-^h* + i^h*-... 



Now add the two series to get 



.2^ , o/ (4) W,4 , / (6) (^),6 



f(x + h) + f(x-h) = 2f(x) + h z f"{x) + 2 J -^h' i + 2 J —^h!> + ... 



Then let 



f( x + h)-2f(x) + f(x-h) t „. . , J^(x) u2 , jW(x) u4 



m = - — ^2 " = f"{x) + r-^h 2 + 2 J -^h A + ..., 



oo 

2fc 



fc=i 



Thus we can use Richardson Extrapolation on ip(h) to get higher order approximations. 
This derivation also gives us the centered difference approximation to the second derivative: 



nx)= f(* + h)-2f(x) + f(x- h ) +o{h2y 



(7.3) 



7.2 Richardson Extrapolation 



The centered difference approximation gives a truncation error of O (/i 2 ) , which is better than 
O (h) . Can we do better? Let's define 

m = ^[f(x+h)-f(x-h)}. 

Had we expanded the Taylor's Series for f(x + h), f(x — h) to more terms we would have seen 
that 

0(h) = f{x) + a 2 h 2 + a 4 h 4 + a 6 h 6 + a s h s + ... 

The constants ctj are a function of f^ +1 \x) only. (In fact, they should take the value of ■) 
What happens if we now calculate 4>(h/2)7 

<f>(h/2) = f'( x ) + ^a 2 h 2 + ^a 4 h 4 + ^a 6 h e + ^a s h 8 + ... 

But we can combine this with (f>{h) to get better accuracy. We have to be a little tricky, but we 
can get the O (/i 2 ) terms to cancel by taking the right multiples of these two approximations: 

4,(h) - = -3f(x) + 1 4h* + -a fj k l> + -a a h* + ... 

3 v 1 4 16 64 
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This approximation has a truncation error of O (/i 4 ) . 

This technique of getting better approximations is known as the Richardson Extrapolation, and 
can be repeatedly applied. We will also use this technique later to get better quadrature rules-that 
is, ways of approximating the definite integral of a function. 



7.2.1 Abstracting Richardson's Method 

We now discuss Richardson's Method in a more abstract framework. Suppose you want to calculate 
some quantity L, and have found, through theory, some approximation: 



<j>{h) = L + a 2kh 



2k 



k=l 



Let 



Now define 



D(n,0) 



£(n, m) 



4 m £(n, m - 1) - £(n - 1, m - 1) 
4m _ J 



(7.4) 



We will be interested in calculating D(n,n) for some n. We claim that 

D{n,n) = L + o(h 2(n+ Vy 

First we examine the recurrence for £(n, m). As in divided differences, we use a pyramid table: 

0(0,0) 

0(1,0) 0(1,1) 

£(2,0) £(2,1) £(2,2) 



£(n,0) £(n, 1) £(ra,2) 



£(n, n) 



By definition we know how to calculate the first column of this table; every other entry in the table 
depends on two other entries, one directly to the left, and the other to the left and up one space. 
Thus to calculate £(n,n) we have to compute this whole lower triangular array. 

We want to show that £(n, n) = L + O (/i 2 ( n+1 )) , that is £(n, n) is a O (/i 2(n+1 )) approximation 
to L. The following theorem gives this result: 

Theorem 7.3 (Richardson Extrapolation). There are constants ak, m such that 

00 / h\ 2k 

D(n,m) = L+ ^ a k ^ m ( — J (0 < m < n) . 

k=m+l ^ ' 



The proof is by an easy, but tedious, induction. We skip the proof. 
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7.2.2 Using Richardson Extrapolation 

We now try out the technique on an example or two. 



Example Problem 7.4. Approximate the derivative of f(x) = logx at x = 1. 

Solution: The real answer is /'(l) = 1/1 = 1, but our computer doesn't know that. Define 



*W-^[«i + *)-/<i-*)l = ^. 

Let's use h = 0.1. We now try to find D(2, 2), which is supposed to be a O (/i 6 ) approximation to 
/'(I) = 1: 



n\m 





1 


2 





« 1.003353477 






1 


log ^l « 1.000834586 


« 0.999994954 




2 


g l « 75 « 1.000208411 


« 0.999999686 


» 1.000000002 



This shows that the Richardson method is pretty good. However, notice that for this simple 
example, we have, already, that 0(0.00001) w 0.999999999. H 

Example Problem 7.5. Consider the ugly function: 

f{x) = arctan(x). 

Attempt to find f'{y/2). Recall that fix) = j^i, so the value that we are seeking is |. 
Solution: Let's use h = 0.01. We now try to find D(2, 2), which is supposed to be a O (/i 6 ) 
approximation to |: 



n\m 





1 


2 





0.333339506181068 






1 


0.333334876543723 


0.333333333331274 




2 


0.33333371913582 


0.333333333333186 


0.333333333333313 



Note that we have some motivation to use Richardson's method in this case: If we let 



m = ^[f(V2+h)-f(V2-h) 



then making h small gives a good approximation to /' (V2~) until subtractive cancelling takes over. 
The following table illustrates this: 
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h 


<Kh) 


1.0 




0.39269908169872408 


0.1 




0.33395069677431943 


0.01 




0.33333950618106845 


0.001 




0.33333339506169679 


0.0001 


0.33333333395058062 


1 x 10" 


-5 


0.33333333334106813 


1 x 10" 


-6 


0.33333333332441484 


1 x 10" 


-7 


0.33333333315788138 


1 x 10" 


-8 


0.33333332760676626 


1 x 10" 


-9 


0.33333336091345694 


1 x 10" 


10 


0.333333360913457 


1 x lO - 


11 


0.333333360913457 


1 x 10- 


12 


0.33339997429493451 


1 x 10- 


13 


0.33306690738754696 


1 x 10- 


14 


0.33306690738754696 


1 x 10- 


15 


0.33306690738754691 


1 x 10- 


16 






The data are illustrated in Figure 7.1. Notice that <j)(h) gives at most 10 decimal places of 
accuracy, then begins to deteriorate; Note however, we get 13 decimal places from D(2,2). H 




1e-12 1 ■ 1 ■ ' ■ ' ■ ' ■ ' ■ ' ■ ' ■ 1 

1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1 



Figure 7.1: The total error for the centered difference approximation to f'(y/2) is shown versus 
h. The total error is the sum of a truncation term which decreases as h decreases, and a roundoff 
term which increases. The optimal h value is around 1 x 10 -5 . Note that Richardson's D(2, 2) 
approximation with h = 0.01 gives much better results than this optimal h. 
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Exercises 



(7.1) Derive the approximation 



/'(*) 



4/ (x + h)- 3/(x) - fjx - 2h) 
6h 



(7.2) 



using Taylor's Theorem. 

(a) Assuming that f(x) has bounded derivatives, give the accuracy of the above approxi- 
mation. Your answer should be something like O (h 7 ). 

(b) Let f(x) = x 3 . Approximate /'(0) with this approximation, using h = \. 

Let f(x) be an analytic function, i.e., one which is infinitely differentiable. Let ip(h) be the 
centered difference approximation to the first derivative: 



f( x + h)-f(x-h) 
2h 



(a) Show that 4>(h) = f'(x) + f f"'(x) + ^f^(x) + ^f^(x) + . . 

(b) Show that 

8 (« W -»/2|) =rw+0(jI)i 



(7.3) Derive the approximation 



h 2 



4/(x + 3h) + 5f(x) - 9fjx - 2h) 
30h 



using Taylor's Theorem. 

(a) What order approximation is this? (Assume f(x) has bounded derivatives of arbitrary 
order.) 

(b) Use this formula to approximate /'(0), where f(x) = x 4 , and h = 0.1 

(7.4) Suppose you want to know quantity Q, and can approximate it with some formula, say (f)(h), 
which depends on parameter h, and such that <p(h) = Q + a±h + a2h 2 + a^h 3 + a^h 4 + . . . 
Find some linear combination of (f>(h) and <j>{—h) which is a O (/i 2 ) approximation to Q. 

(7.5) Assuming that 4>(h) = Q + a2h 2 + a^h 4 + aeh 6 . . ., find some combination of 4>(h), 4>(h/3) 
which is a O (/i 4 ) approximation to Q. 

(7.6) Let A be some number in (0, 1). Assuming that <j)(h) = Q + a^l? + a^h 4 + agh 6 . . ., find some 
combination of (f)(h) , 4>(Xh) which is a O (h 4 ) approximation to Q. To make the constant 
associated with the h 4 term small in magnitude, what should you do with A? Is this practical? 
Note that the method of Richardson Extrapolation that we considered used the value A = 1/2. 

(7.7) Assuming that <p{h) = Q + a2h 2 + a^h 4 + a^h 6 . . ., find some combination of (p(h),(j>(h/4) 
which is a O (h 4 ) approximation to Q. 

(7.8) Suppose you have some great computational approximation to the quantity Q such that 
ij)(h) = Q + a%h 3 + a^h 6 + agh 9 . . . Can you find some combination of ip(h),ip(h/2) which is 
a O (/i 6 ) approximation to Ql 

(7.9) Complete the following Richardson's Extrapolation Table, assuming the first column consists 
of values D(n, 0) for n = 0, 1, 2: 



n\m 





1 


2 





2 






1 


1.5 


? 




2 


1.25 


? 


? 



(See equation 7.4 if you've forgotten the definitions.) 
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10) Write code to complete a Richardson's Method table, given the first column. Your m-file 
should have header line like: 
function Dnn = richardsons (colO) 

where Dnn is the value at the lower left corner of the table, D(n, n) while colO is the column 
of n + 1 values D(i, 0), for i = 0, 1, . . . , n. Test your code on the following input: 
octave :1> colO = [1 0.5 0.25 0.125 0.0625 0.03125]; 
octave :2> richardsons (colO) 
ans = 0.019042 

(a) What do you get when you try the following? 

octave :5> colO = [1.5 0.5 1.5 0.5 1.5 0.5 1.5]; 
octave :6> richardsons (colO) 

(b) What do you get when you try the following? 

octave :7> colO = [0.9 0.99 0.999 0.9999 0.99999]; 
octave :8> richardsons (colO) 



Chapter 8 

Integrals and Quadrature 



8.1 The Definite Integral 



Often enough the numerical analyst is presented with the challenge of finding the definite integral 
of some function: 



In your golden years of Calculus, you learned the Fundamental Theorem of Calculus, which claims 
that if /(x) is continuous, and F(x) is an antiderivative of f(x), then 



What you might not have been told in Calculus is there are some functions for which a closed 
form antiderivative does not exist or at least is not known to humankind. Nevertheless, you may 
find yourself in a situation where you have to evaluate an integral for just such an integrand. An 
approximation will have to do. 

8.1.1 Upper and Lower Sums 

We will review the definition of the Riemann integral of a function. A partition of an interval [a, b] 
is a finite, ordered collection of nodes xf. 



Given such a partition, P, define the upper and lower bounds on each subinterval [ 
follows: 





a = xq < x\ < X2 < ■ ■ ■ < x. 



n 



= b. 



mi = inf {/(x) | Xi < x < x i+1 } 
Mi = sup {/(x) | Xi < x < x i+1 } 



Then for this function / and partition P, define the upper and lower sums: 



n-1 




i=0 



n-1 



U(f,P) 



^2 Mi (x i+ i - x^ 



i=0 
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We can interpret the upper and lower sums graphically as the sums of areas of rectangles defined 
by the function / and the partition P, as in Figure 8.1. 




(a) The Lower Sum 



(b) The Upper Sum 



Figure 8.1: The (a) lower, and (b) upper sums of a function on a given interval are shown. These 
approximations to the integral are the sums of areas of rectangles. Note that the lower sums are 
an underestimate, and the upper sums an overestimate of the integral. 

Notice a few things about the upper, lower sums: 



(ii) If we switch to a "better" partition (i.e., a finer one), we expect that L(/, •) increases and 
U(f,-) decreases. 

The notion of integrability familiar from Calculus class (that is Riemann Integrability) is defined 
in terms of the upper and lower sums. 

Definition 8.1. A function / is Riemann Integrable over interval [a, b] if 



where the supremum and infimum are over all partitions of the interval [a, b\. Moreover, in case 
f(x) is integrable, we define the integral 



You may recall the following 

Theorem 8.2. Every continuous function on a closed bounded interval of the real line is Riemann 
Integrable (on that interval). 

Continuity is sufficient, but not necessary. 

Example 8.3. Consider the Heaviside function: 



This function is not continuous on any interval containing 0, but is Riemann Integrable on every 
closed bounded interval. 



(i) L(f,P)<U(f,P). 



supL(/,P)=inf£/(/,P) 
p P 



P 
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Example 8.4. Consider the Dirichlet function: 

. . _ f x rational 
} 1 i irrational 

For any partition P of any interval [a, b], we have L(f, P) = 0, while U (/, P) = 1, so 

supL(/,P) = 0^1 = inf£/(/,P), 
p p 

so this function is not Riemann Integrable. 
8.1.2 Approximating the Integral 

The definition of the integral gives a simple method of approximating an integral f(x) dx. The 
method cuts the interval [a, b] into a partition of n equal subintervals X{ = a+^p, for i = 0, 1, . . . , n. 
The algorithm then has to somehow find the supremum and infimum of f{x) on each interval 
[xj,Xj + i]. The integral is then approximated by the mean of the lower and upper sums: 

cb 



J f{x)6xKt\{L{f,P) + U{f,P)) 



Because the value of the integral is between L(f,P) and U(f,P), this approximation has error at 
most 

\ (U(f, P)~L(f,P)). 

Note that in general, or for a black box function, it is usually not feasible to find the suprema 
and infima of f(x) on the subintervals, and thus the lower and upper sums cannot be calculated. 
However, if some information is known about the function, it becomes easier: 

Example 8.5. Consider for example, using this method on some function f(x) which is monotone 
increasing, that is x < y implies f(x) < f(y). In this case, the infimum of f(x) on each interval 
occurs at the leftmost endpoint, while the supremum occurs at the right hand endpoint. Thus for 
this partition, P, we have 

n— 1 I n— 1 

£(/> P ) = ^2 m i \ x k+l ~ Xk\ = ^2 f( X k) 

k=0 k=0 

n—l |, | n— 1 i, I n 

u(f, p) = J2^ - x k\ = E f( x ^ = E 

k=0 k=0 k=l 

Then the error of the approximation is 
8.1.3 Simple and Composite Rules 

For the remainder of this chapter we will study "simple" quadrature rules, i.e., rules which approx- 
imate the integral of a function, f(x) over an interval [a, b] by means of a number of evaluations of 
/ at points in this interval. The error of a simple quadrature rule usually depends on the function 
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f, and the width of the interval [a, b] to some power which is determined by the rule. That is we 
usually think of a simple rule as being applied to a small interval. 

To use a simple rule on a larger interval, we usually cast it into a "composite" rule. Thus the 
trapezoidal rule, which we will study next becomes the composite trapezoidal rule. The means of 
extending a simple rule to a composite rule is straightforward: Partition the given interval into 
subintervals, apply the simple rule to each subinterval, and sum the results. Thus, for example if 
the interval in question is [a, 0\, and the partition is a = xq < x\ < X2 < ■ ■ ■ < x n = j3, we have 

n-1 

composite rule on [a, (3] = simple rule applied to 

i=0 

8.2 Trapezoidal Rule 

Suppose we are trying to approximate the integral 




f{x) dx, 



for some unpleasant or black box function f(x). 
The trapezoidal rule approximates the integral 

/ f{x) dx 

J a 

by the (signed) area of the trapezoid through the points (a, /(a)) , (6, f(b)) , and with one side the 
segment from a to b. See Figure 8.2. 




Figure 8.2: The trapezoidal rule for approximating the integral of a function over [a, b] is shown. 
By old school math, we can find this signed area easily. This gives the (simple) trapezoidal rule: 



8.2. TRAPEZOIDAL RULE 
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The composite trapezoidal rule can be written in a simplified form, one which you saw in your 
calculus class, if the interval in question is partitioned into equal width subintervals. That is if you 
let [a, P] be partitioned by 

a = xq < x\ < X2 < x 3 < . . . < x n = (3, 
with Xi = a + ih, where h = ((3 — a)/n, then the composite trapezoidal rule is 

/ f{x) dx = £ / f(?) ~ ? E (^+1 " x *) \f( x *) + • t 8 - 1 ) 

1/0 i=0 ^ i=0 

Since each subinterval has equal width, Xi + \ — x% = h, and we have 



(8.2) 




In your calculus class, you saw this in the less comprehensible form: 



["fix) 

J a 



dx « h 



f(a) + f(b) 



n-l 



i=i 



Note that the composite trapezoidal rule for equal subintervals is the same as the approximation 
we found for increasing functions in Example 8.5. 



Example Problem 8.6. Approximate the integral 



f— 

Jo 1+: 



dx 



by the composite trapezoidal rule with a partition of equally spaced points, for n = 2. 

Solution: We have h = = 1, and /(xo) = l,/(xi) = ^,/(x2) = \. Then the composite 

trapezoidal rule gives the value 



I[ / (x ) + /(x 1 ) + /(x 1 ) + /(x 2 )] = i 



1 + 1 + 



1 



11 

10' 



The actual value is arctan2 « 1.107149, and our approximation is correct to two decimal places. H 



8.2.1 How Good is the Composite Trapezoidal Rule? 

We consider the composite trapezoidal rule for partitions of equal subintervals. Let Pi(x) be the 
polynomial of degree < 1 that interpolates f(x) at Xj,Xj-|-i. Let 

/i = /(*) dx, Tj = f i+1 Pi (x) dx = (xi+i - x.) ±^!±l) = £ {f{xi) + 

That's right: the composite trapezoidal rule approximates the integral of fix) over [xj,Xj+i] by 
the integral of pi(x) over the same interval. 

Now recall our theorem on polynomial interpolation error. For x € [xj,Xj+i] , we have 

f(x) - Pi{x) = — /( 2 )(£,.) (x - ^i) (a: - x i+1 ) , 
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for some £ x G [x«, Xj+i] . Recall that ^ depends on x. To make things simpler, call it 
Now integrate: 

h-Ti= I f(x) - pi(x) dx = - / f"(€(x)) (x - Xi) (x - x i+1 ) dx. 

J Xi ^ J Xi 

We will now attack the integral on the right hand side. Recall the following theorem: 

Theorem 8.7 (Mean Value Theorem for Integrals). Suppose / is continuous, g is Riemann 
Integrable and does not change sign on [a, 0\. Then there is some ( G [a,0\ such that 

f(x)g(x) dx = /(C) / g(x) dx. 



We use this theorem on our integral. Note that (x — Xi) (x — Xj+i) is nonpositive on the interval 
of question, [xj,Xj+i]. We assume continuity of f"(x), and wave our hands to get continuity of 
f"(£(x)). Then we have 



h-Ti = \f"{£) [ (x - Xi) {x - x i+1 ) dx, 

Z J Xi 



for some £j G By boring calculus and algebra, we find that 



JXi 



Xi+i /j3 

(x - Xi) (x - x i+ i) dx = - — . 



This gives 

h 3 

I'i ~ Ti = ~~^f 

for some & G [xi,x i+1 ]. 

We now sum over all subintervals to find the total error of the composite trapezoidal rule 



n— 1 , 3 ra— 1 /, \ 7 2 



12 K ' 12 

i=0 i=0 



1 71— 1 

-E/"(6 



i=0 

On the far right we have an average value, ^ Y17=o /"(£»)) which lies between the least and greatest 
values of /" on the inteval [a, b], and thus by the IVT, there is some £ which takes this value. So 



This gives us the theorem: 



Theorem 8.8 (Error of the Composite Trapezoidal Rule). Let f"(x) be continuous on [a, b]. 
Let T be the value of the trapezoidal rule applied to f(x) on this interval with a partition of uniform 
spacing, h, and let / = f(x) dx. Then there is some £ G [a, b] such that 

Note that this theorem tells us not only the magnitude of the error, but the sign as well. Thus 
if, for example, f(x) is concave up and thus /" is positive, then I — T will be negative, i.e., the 
trapezoidal rule gives an overestimate of the integral /. See Figure 8.3. 
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Figure 8.3: The trapezoidal rule is an overestimate for a function which is concave up, i.e., has 
positive second derivative. 

8.2.2 Using the Error Bound 

Example Problem 8.9. How many intervals are required to approximate the integral 

In 2 = I = / dx 

Jo 1 + x 

to within 1 x 10~ 10 ? 

Solution: We have f(x) = thus f{x) = -j^- And f"{x) = Thus /"(£) is 

continuous and bounded by 2 on [0, 1]. If we use n equal subintervals then Theorem 8.8 tells us the 
error will be 

i-o/i-o\ 2 no 

12 \ n J 1 ^' 12n 2 ' 
To make this smaller than 1 x 10~ 10 , in absolute value, we need only take 

i<lx 10- 10 , 

and so n > x 10 5 suffices. Because f"{x) is positive on this interval, the trapezoidal rule will 
be an overestimate. H 

Example Problem 8.10. How many intervals are required to approximate the integral 

x 3 - 1 dx 



o 

to within 1 x 10 -6 ? 

Solution: We have f(x) = x 3 — 1, thus f'{x) = 3x 2 , and f"(x) = 6x. Thus f"(£) is continuous 
and bounded by 12 on [0, 2]. If we use n equal subintervals then by Theorem 8.8 the error will be 

2-0 (2-Q\ 2 „ 2f"(0 



12 V n ) 3n 2 
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To make this smaller than 1 x 10 6 , in absolute value, it suffices to take 

24 « 

s? 5 1 x 10 

and so n > \/8 x 10 3 suffices. Because f"(x) is positive on this interval, the trapezoidal rule will 
be an overestimate. H 



8.3 Romberg Algorithm 

Theorem 8.8 tells us, approximately, that the error of the composite trapezoidal rule approximation 
is O (h 2 ). If we halve h, the error is quartered. Sometimes we want to do better than this. We'll 
use the same trick that we did from Richardson extrapolation. In fact, the forms are exactly the 
same. 

Towards this end, suppose that /, a, b are given. For a given n, we are going to use the trapezoidal 
rule on a partition of 2 n equal subintervals of [a, b}. That is h = ^r- Then define 

lb - a 2 ™ _1 

i=0 



a 



2 n 



f(a) f(b) 2 ^\f .b-a 



The intervals used to calculate 4>(n + 1) are half the size of those for (f)(n). As mentioned above, 
this means the error is one quarter. 

It turns out that if we had proved the error theorem differently, we would have proved the 
relation: 

rb 

<K n ) = / f( x ) dx + a 2 h 2 l + a 4 /4 + a^h® + a%h % n + . . . , 

J a 

where h n = ^pr- The constants are a function of f^\x) only. This should look just like something 
from Chapter 7. What happens if we now calculate 4>(n + 1)? We have 

rb 

(f)(n + 1) = / f(x) dx + a 2 h 2 n+1 + a A h A n+1 + a 6 /i^ +1 + a 8 ^ +1 + . . . , 

J a 

rb i 11 1 

= J f(x) dx + -a 2 h 2 n + — a 4 h n + -^a G h % n + ^a 8 /i® + . . . . 

This happens because h n+ \ = = = \- As with Richardon's method for approximating 

derivatives, we now combine the right multiples of these: 



/ 3 15 63 

fix) dx + -4/4 + — a % h & n + — a 8 h 8 n + ... 

/b 15 21 

f{x) dx - -4/4 - r^a 6 h b n - ^a 8 hf 



4o(n)-0(n+\) _ , , (i , ,,. , 



This approximation has a truncation error of O (/i^) . 

Like in Richardson's method, we can use this to get better and better approximations to the 
integral. We do this by constructing a triangular array of approximations, each entry depending 
on two others. Towards this end, we let 

R(n,0) = (f)(n), 
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then define, for m > 

A m F?(n rn — 1\ — Ff(n — 1 m — 1\ 

(8.3) 



, 4 m R(n, m - 1) - R(n - 1, m - 1) 
R(n,m) = 



4m _ 1 



The familiar pyramid table then is: 

12(0,0) 

i?(l,0) i?(l,l) 

R(2,0) R(2,l) R(2,2) 

R(n,0) R(n,l) R(n,2) ■■■ R(n,n) 

Even though this is exactly the same as Richardon's method, it has another name: this is called 
the Romberg Algorithm. 

Example Problem 8.11. Approximating the integral 

dx 



o 1 + ^ 

by Romberg's Algorithm; find R(l, 1). 

Solution: The first column is calculated by the trapezoidal rule. Successive columns are found by 
combining members of previous columns. So we first calculate R(0, 0) and R(l, 0). These are fairly 
simple, the first is the trapezoidal rule on a single subinterval, the second is the trapezoidal rule on 
two subintervals. Then 

fl(0,0) = ^^[/(0) + /(2)] = jj, 

R(l, 0) = 2 -^\ [/(0) + /(l) + /(l) + /(2)] = 

Then, using Romberg's Algorithm we have 

4fl(l,0)-i i (0,0) = t-j| = 32 = 
v ' ; 4-1 3 30 

H 

At this point we are tempted to use Richardson's analysis. This would claim that R(n, n) is a 
O (h^ n+1 ^ approximation to the integral. However, ho = b — a, and need not be smaller than 1. 

This is a bit different from Richardson's method, where the original h is independently set before 
starting the triangular array; for Romberg's algorithm, ho is determined by a and b. 

We can easily deal with this problem by picking some k such that is small enough, say 
smaller than 1. Then calculating the following array: 

R(k,0) 
R(k + 1,0) R{k + 1,1) 
R{k + 2,0) J?(ife + 2,1) R(k + 2,2) 



R(k + n,0) R(k + n,l) R(k + n,2) ••• R(k + n,n) 
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Quite often Romberg's Algorithm is used to compute columns of this array. Subtractive can- 
celling or unbounded higher derivatives of f(x) can make successive approximations less accurate. 
For this reason, entries in ever rightward columns are usually not calculated, rather lower entries 
in a single column are calculated instead. That is, the user calculates the array: 



R(k,0) 
R(k + 1,0) 
R(k + 2,0) 



R(k + 1,1) 
R(k + 2,1) 



R(k + 2,2) 



R(k + n,0) R(k + n,l) R(k + n,2) 

R(k + n + 1,0) R(k + n + 1,1) R{k + n + 1,2) 

R(k + n + 2,0) R(k + n + 2,1) R(k + n + 2,2) 

R{k + n + 3,0) R(k + n + 3,l) R(k + n + 3,2) 



R(k + n, n) 
R{k + n + l,n) 
R(k + n + 2,n) 
R(k + n + 3,n) 



Then R(k + n + l,n) makes a fine approximation to the integral as I — > oo. Usually n is small, 
like 2 or 3. 



8.3.1 Recursive Trapezoidal Rule 

It turns out there is an efficient way of calculating R(n + 1, 0) given R(n, 0); first notice from the 
above example that 

R(Q,0) = t^Lh f{a) + mh 



R(i,o) = 



1 2 
b-a l 

2 2 



f(a) + f(^) + f(^) + f(b) 



It would be best to calculate R(l, 0) without recalculating f(a) and f{b). It turns out this is possible. 



Let h n = and recall that 



R{n, 0) = 4>{n) = h n 

Thus 

R(n + 1,0) = <t>(n + 1) = h n+1 



/(«) + /(&) 



2 n -l 



+ Yl f( a + ih n) 



i=l 



2 n + l_ 1 



1. 

= -h, 



2 + Z> f{a + ih n+1 ) 

/(«) + f(b) 



2 n -l 



+ J2f(a + (2i- l)h n+1 ) + / (a + (2i)h n ^j 



/(«) + f(b) 



2 n -l 



2 n -l 



+ f (a + */*n) + / (a + (2i - l)/i n+1 ) 



2 n -l 



]-R{n, 0) + /i n+ i J] / (a + (2i - l)h n+l ) 



i=i 



Then calculating + 1,0) requires only 2™ — 1 additional evaluations of f{x), instead of the 
2 n+l + 1 usually required. 
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8.4 Gaussian Quadrature 

The word quadrature refers to a method of approximating the integral of a function as the linear 
combination of the function at certain points, i.e., 



/ f{x) dx « A f(x ) + A 1 f(x 1 ) + ... A n f{x n ), 

J a 

for some collection of nodes {xj}™ =0 , and weights {Ai}™ =Q . Normally one finds the nodes and weights 
in a table somewhere; we expect a quadrature rule with more nodes to be more accurate in some 
sense-the tradeoff is in the number of evaluations of /(■)■ We will examine how these rules are 
created. 



(8.4) 



8.4.1 Determining Weights (Lagrange Polynomial Method) 

Suppose that the nodes {xj}™ =0 are given. An easy way to find "good" weights {^4j}™ =0 for these 
nodes is to rig them so the quadrature rule gives the integral of p(x), the polynomial of degree < n 
which interpolates f[x) on these nodes. Recall 



i=0 

where £i(x) is the i th Lagrange polynomial. Thus our rigged approximation is the one that gives 

rb rb n rb 

j f(x)dx& / p(x) dx = y^J(xj) / li{x)dx. 



i=0 

If we let 

r b 

A. 



= / £i{x)dx, 

J a 



then we have a quadrature rule. 

If f(x) is a polynomial of degree < n then f(x) = p(x), and the quadrature rule is exact. 

Example Problem 8.12. Construct a quadrature rule on the interval [0,4] using nodes 0, 1, 2. 
Solution: The nodes are given, we determine the weights by constructing the Lagrange Polyno- 
mials, and integrating them. 

4>(S) = ^~^n"o' = ^ X ~ " 2 )' 

'-..(-) - "rzr:: --ix)i,-2). 



(x 


-l)(x 


-2) 


(0 


-1)(0 


-2) 


(x 


- 0)(x 


-2) 


(1 


-0)(1 


-2) 


(x 


- 0)(x 


-1) 


(2 


-0)(2 


-1) 



^(x) = ) n n ;; o ^ = -(x)(x - 1). 



Then the weights are 

A) = J to(x) dx = J i(x - l)(x - 2) dx = |, 

f 4 f A 16 

A 1 = / £ 1 (x)dx= / -(x)(x-2)dx = , 

Jo Jo 3 

A 2 = £2(x)dx= / -(x)(x-l)dx 

Jo Jo 2 



20 
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Thus our quadrature rule is 



f 8 16 20 

^ /(x)dx«-/(0)- T /(l) + T /(2). 



We expect this rule to be exact for a quadratic function f(x). To illustrate this, let f(x) 
By calculus we have 



x 2 + l. 



Jo 



x 2 + 1 dx= -x 3 + x 



64 „ 76 

y +4 = y 



The approximation is 



£, 2 + ld*^[0 + l]-^[l + l] + f [4 + 1] = f. 



8.4.2 Determining Weights (Method of Undetermined Coefficients) 

Using the Lagrange Polynomial Method to find the weights Ai is fine for a computer, but can 
be tedious (and error-prone) if done by hand (say, on an exam). The method of undetermined 
coefficients is a good alternative for finding the weights by hand, and for small n. 

The idea behind the method is to find n + 1 equations involving the n + 1 weights. The equations 
are derived by letting the quadrature rule be exact for f(x) = x J for j = 0, 1, . . . , n. That is, setting 

rb n 

/ x j dx = ^ AiixiY . 
Ja i=0 

For example, we reconsider Example Problem 8.12. 

Example Problem 8.13. Construct a quadrature rule on the interval [0,4] using nodes 0, 1, 2. 
Solution: The method of undetermined coefficients gives the equations: 

J 1 dx = 4 = A + A 1 + A 2 
Jo 

4 

x dx = 8 = A-y + 2A 2 

4 

x 2 dx = 64/3 = Ai +AA 2 . 
We perform Naive Gaussian Elimination on the system: 





1 


1 




(• 


1 


2 






1 


4 


64/3 / 



We get the same weights as in Example Problem 8.12: A 2 = ^-,Ai = — A) = §• ^ 

Notice the difference compared to the Lagrange Polynomial Method: undetermined coefficients 
requires solution of a linear system, while the former method calculates the weights "directly." 
Since we will not consider n to be very large, solving the linear system may not be too burdensome. 

Moreover, the method of undetermined coefficients is useful in more general settings, as illus- 
trated by the next example: 



L 



L 
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Example Problem 8.14. Determine a "quadrature" rule of the form 

f 1 f(x)dx^Af(l) + Bf'(l) + Cf"(l) 
Jo 

that is exact for polynomials of highest possible degree. What is the highest degree polynomial for 
which this rule is exact? 

Solution: Since there are three unknown coefficients to be determined, we look for three equations. 
We get these equations by plugging in successive polynomials. That is, we plug in f(x) = 1, x, x 2 
and assuming the coefficients give equality: 

f Idx = 1 = Al +B0+C0 = A 
Jo 

/ xdx = l/2 = Al + Bl + CO = A + B 
Jo 

[ x 2 dx = l/3 = Al + B2 + C2 = A + 2B + 2C 
Jo 

This is solved by A = 1, B = — 1/2, C = 1/6. This rule should be exact for polynomials of degree 
no greater than 2, but it might be better. We should check: 

\ x 3 dx = 1/4/ 1/2 = 1-3/2 + 1 = A1 + B3 + C6, 
Jo 

and thus the rule is not exact for cubic polynomials, or those of higher degree. H 



8.4.3 Gaussian Nodes 

It would seem this is the best we can do: using n + 1 nodes we can devise a quadrature rule that is 
exact for polynomials of degree < n by choosing the weights correctly. It turns out that by choosing 
the nodes in the right way, we can do far better. Gauss discovered that the right nodes to choose 
are the n + 1 roots of the (nontrivial) polynomial, q(x), of degree n + 1 which has the property 

f b 

/ x k q(x)dx = (0<k<n). 

J a 

(If you view the integral as an inner product, you could say that q(x) is orthogonal to the polyno- 
mials x k in the resultant inner product space, but that's just fancy talk.) 

Suppose that we have such a q(x)-we will not prove existence or uniqueness. Let f(x) be a 
polynomial of degree < 2n + 1. We write 



f(x) =p(x)q(x) +r(x). 
Both p(x),r(x) are of degree < n. Because of how we picked q(x) we have 



J dx = 0. 

J a 

Thus 



/ p(x)q(x)< 

J a 



f'O f'O f'O f'O 

/ f(x) dx = p(x)q(x) dx + / r(x)dx = / r(x)dx. 

J a J a J a J a 
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Now suppose that the Ai are chosen by Lagrange Polynomials so the quadrature rule on the 
nodes Xi is exact for polynomials of degree < n. Then 

n n n 

^2 A if( X i) = ^2 A i [p{Xi)q(xi] + r(Xi)] = ^ A i r ( x i)- 
i=0 i=0 i=0 

The last equality holds because the x\ are the roots of q(x). Because of how the Ai are chosen we 
then have 

n n r-b f-b 

S ^A i f{xi) = ^A i r{xi)=j r(x)dx = f(x)dx. 
i=0 i=0 Ja Ja 

Thus this rule is exact for f(x). We have (or rather, Gauss has) made quadrature twice as good. 

Theorem 8.15 (Gaussian Quadrature Theorem). Let xi be the n + 1 roots of a (nontrivial) 
polynomial, q(x), of degree n + 1 which has the property 

rb 

I x k q{x) dx = (0 < k < n) . 

J a 

Let Ai be the coefficients for these nodes chosen by integrating the Lagrange Polynomials. Then the 
quadrature rule for this choice of nodes and coefficients is exact for polynomials of degree < 2n + 1. 

8.4.4 Determining Gaussian Nodes 

We can determine the Gaussian nodes in the same way we determine coefficients. The example is 
illustrative 

Example Problem 8.16. Find the two Gaussian nodes for a quadrature rule on the interval [0, 2]. 
Solution: We will find the function q(x) of degree 2, which is "orthogonal" to l,x under the inner 
product of integration over [0,2]. Thus we let q(x) = cq + c\x + C2X 2 . The orthogonality condition 
becomes 

,2 ,2 

/ lq(x) dx = xq(x) dx = that is, 

Jo Jo 

/ Co + c\x + C2X 2 dx = cqx + c\x 2 + C2X 3 dx = 
Jo Jo 

Evaluating these integrals gives the following system of linear equations: 

8 

2c + 2ci + -c 2 = 0, 
g 

2c + -ci + 4c 2 = 0. 

This system is "under determined," that is, there are two equations, but three unknowns. Notice, 
however, that if q(x) satisfies the orthogonality conditions, then so does q(x) = aq(x), for any real 
number a. That is, we can pick the scaling of q(x) as we wish. 

With great foresight, we "guess" that we want c 2 = 3. This reduces the equations to 

2c + 2ci = -8, 
2c + -ci = -12. 
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Simple Gaussian Elimination (cf. Chapter 3) yields the answer cq = 2, c\ = —6, C2 = 3. 
Then our nodes are the roots of q(x) = 2 — Qx + 3x 2 . That is the roots 

6 ± V36 - 24 ^3 
6 ~ ~3~' 

These nodes are a bit ugly. Rather than construct the Lagrange Polynomials, we will use the 
method of undetermined coefficients. Remember, we want to construct Aq,A\ such that 



L 



2 f(x) dx « A f(l - ^) + A 1 f(l + ^) 



is exact for polynomial /(x) of degree < 1. It suffices to make this approximation exact for the 
"building blocks" of such polynomials, that is, for the functions 1 and x. That is, it suffices to find 
A ,Ai such that 



s: 

L 



1 dx = A + A 1 
2 xdx = A (l-^-) + A 1 (l + ^-) 



This gives the equations 



2 = A + A x 

2 = A)(l-^) + A 1 (l + ^) 



This is solved by A§ = A\ = 1. 
Thus our quadrature rule is 



^/(x)dx«/(l-^) + /(l + ^ 



V3, 



We expect this rule to be exact for cubic polynomials. 

Example Problem 8.17. Verify the results of the previous example problem for f(x) 
Solution: We have 



o 



2 2 

f{x)dx = (1/4) x 4 



4. 



The quadrature rule gives 

/(I- f ) + /(! + f) 



3 



+ 1 



3 



1 _j^ + 8 ?-?^U[l + 8^ + 8? + ^ 
3 9 27 3 9 27 



= (2 - V3 - \/3/9) + (2 + + \/3/9) = 4 
Thus the quadrature rule is exact for f(x). 
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8.4.5 Reinventing the Wheel 

While it is good to know the theory, it doesn't make sense in practice to recompute these things all 
the time. There are books full of quadrature rules; any good textbook will list a few. The simplest 
ones are given in Table 8.1. 

See also http : //mathworld . wolfram . com/Legendre-GaussQuadrature . html 



n 




Ai 





x = 


A = 2 


1 


Xq = -y/1/3 


Aq = 1 


Xl = a/173 


A 1 = l 




xq = -a/ 3 / 5 


A = 5/9 


2 


xi = 


A x = 8/9 




x 2 = a/3/5 


Ay = 5/9 



Table 8.1: Gaussian Quadrature rules for the interval [—1,1]. Thus f^ 1 f(x)dx « Yl^o Af{xi), 
with this relation holding exactly for all polynomials of degree no greater than 2n + 1. 

Quadrature rules are normally given for the interval [—1,1]. On first consideration, it would 
seem you need a different rule for each interval [a, b]. This is not the case, as the following example 
problem illustrates: 

Example Problem 8.18. Given a quadrature rule which is good on the interval [—1, 1], derive a 
version of the rule to apply to the interval [a, b\. 
Solution: Consider the substitution: 



Then 



b — a b + a b — a 

x = 1-\ , so ax = at. 

2 2 2 

/*^ . . f 1 „ fb-a b + a\ b-a , 



Letting 



2 2 

if f(x) is a polynomial, g(t) is a polynomial of the same degree. Thus we can use the quadrature 
rule for [—1, 1] on g(t) to evaluate the integral of f(x) over [a, b\. H 

Example Problem 8.19. Derive the quadrature rules of Example Problem 8.16 by using the 
technique of Example Problem 8.18 and the quadrature rules of Table 8.1. 
Solution: We have a = 0, b = 2. Thus 

b — a„(b — a b + a\ „, „. 
rf*) = — /(—*+— )=/(* + !)■ 

To integrate f(x) over [0, 2], we integrate g(t) over [—1, 1]. The standard Gaussian Quadrature rule 
approximates this as 

1 g(t) dt « g(- a/T73) + g(V^) = /(I " + /(I + V 7 ^)- 

-l 

This is the same rule that was derived (with much more work) in Example Problem 8.16. H 
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Exercises 

(8.1) Use the composite trapezoidal rule, by hand, to approximate 

r 3 

/ x 2 dx (= 9) 
Jo 

Use the partition {xi} 2 =0 = {0, 1,3} . Why is your approximation an overestimate? 

(8.2) Use the composite trapezoidal rule, by hand, to approximate 



1 i 



X + 1 



dx (= In 2 » 0.693) 



Use the partition {xi} 3 =0 = {0, \, ^, l} . Why is your approximation an overestimate? (Check: 
I think the answer is 0.7) 
(8.3) Use the composite trapezoidal rule, by hand to approximate 



L 



1 A 

dx. 



l + x 



2 



Use n = 4 subintervals. How good is your answer? 

(8.4) Use Theorem 8.8 to bound the error of the composite trapezoidal rule approximation of 
J 2 x 3 dx with n = 10 intervals. You should find that the approximation is an overestimate. 

(8.5) How many equal subintervals of [0, 1] are required to approximate Jq 1 cos x dx with error 
smaller than 1 x 10~ 6 by the composite trapezoidal rule? (Use Theorem 8.8.) 

(8.6) How many equal subintervals would be required to approximate 

o dx. 

o l + x 2 

to within 0.0001 by the composite trapezoidal rule? (Hint: Use the fact that \f"(x)\ < 8 on 
[0, 1] for f(x) = 4/(1 + x 2 )) 

(8.7) How many equal subintervals of [2, 3] are required to approximate J* 2 3 e x dx with error smaller 
than 1 x 10~ 3 by the composite trapezoidal rule? 

(8.8) Simpson's Rule for quadrature is given as 



/ 

J a 



A x 

f{x) dx « — [f(x ) + 4/0n) + 2f(x 2 ) + 4/(x 3 ) + • • • + 2/(x n _ 2 ) + 4/(x n _i) + f(x n )] , 



where Ax = (b — a)/n, and n is assumed to be even. Show that Simpson's Rule for n = 2 is 
actually given by Romberg's Algorithm as R(l, 1). As such we expect Simpson's Rule to be 
a O (h A ) approximation to the integral. 
(8.9) Find a quadrature rule of the form 



f 1 f(x) dx « Af(0) + Bf (1/2) + Cf(l) 
Jo 



that is exact for polynomials of highest possible degree. What is the highest degree polynomial 
for which this rule is exact? 
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(8.10) Determine a "quadrature" rule of the form 



/ f(x)dx^Af(0) + Bf'(-l) + Cf'(l) 



J-i 



that is exact for polynomials of highest possible degree. What is the highest degree polynomial 
for which this rule is exact? (Since this rule uses derivatives of /, it does not exactly fit our 
definition of a quadrature rule, but it may be applicable in some situations.) 
(8.11) Determine a "quadrature" rule of the form 



that is exact for polynomials of highest possible degree. What is the highest degree polynomial 
for which this rule is exact? 
(8.12) Consider the so-called order n Chebyshev Quadrature rule: 



Find the weighting c n and nodes Xi for the case n = 2 and the case n = 3. For what order 
polynomials are these rules exact? 
(8.13) Find the Gaussian Quadrature rule with 2 nodes for the interval [1,5], i.e., find a rule 



Before you solve the problem, consider the following questions: do you expect the nodes to 
be the endpoints 1 and 5? do you expect the nodes to be arranged symmetrically around the 
midpoint of the interval? 
(8.14) Find the Gaussian Quadrature rule with 3 nodes for the interval [—1, 1], i.e., find a rule 



To find the nodes xq, x\, X2 you will have to find the zeroes of a cubic equation, which could be 
difficult. However, you may use the simplifying assumption that the nodes are symmetrically 
placed in the interval [—1, 1]. 

(8.15) Write code to approximate the integral of a / on [a, b] by the composite trapezoidal rule on 
n equal subintervals. Your m-file should have header line like: 

function iappx = trapezoidal (f , a, b,n) 
You may wish to use the code: 
x = a .+ (b-a) .* (0:n) ./ n; 

If f is defined to work on vectors element-wise, you can probably speed up your computation 
by using 

bigsum = 0.5 * ( f(x(l)) + f(x(n+l)) ) + sum( f(x(2:(n))) ); 

(8.16) Write code to implement the Gaussian Quadrature rule for n = 2 to integrate / on the 
interval [a, b] . Your m-file should have header line like: 

function iappx = gauss2 (f , a,b) 

(8.17) Write code to implement composite Gaussian Quadrature based on code from the previous 
problem. Something like the following probably works: 
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function iappx = gaussComp(f ,a,b,n) 

% code to approximate integral of f over n equal subintervals of [a,b] 
x = a .+ (b-a) .* (0:n) ./ n; 
iappx = 0; 
for i=l:n 

iappx += gauss2(f ,x(i) ,x(i+l)) ; 
end 

Use your code to approximate the error function: 



Compare your results with the octave/Matlab builtin function erf. (Try help erf in octave, 
or see http : / /mathworld . wolfram . com/Erf . html) 

The error function is used in probability. In particular, the probability that a normal random 
variable is within z standard deviations from its mean is 



Thus erf(l/\/2) ~ 0.683, and erf(2/\/2) ~ 0.955. These numbers should look familiar to you. 




ed(z/y/2) 
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Chapter 9 

Least Squares 



9.1 Least Squares 

Least squares is a general class of methods for fitting observed data to a theoretical model function. 
In the general setting we are given a set of data 



X 


x 


x 1 




X n 


y 


Vo 


yi 




y n 



and some class of functions, T. The goal then is to find the "best" / € T to fit the data to y = f(x). 
Usually the class of functions T will be determined by some small number of parameters; the number 
of parameters will be smaller (usually much smaller) than the number of data points. The theory 
here will be concerned with defining "best," and examining methods for finding the "best" function 
to fit the data. 

9.1.1 The Definition of Ordinary Least Squares 

First we consider the data to be two vectors of length n + 1. That is we let 



x = [xq X\ ... 



and y = [y y 1 ... y n ] T . 



The error, or residual, of a given function with respect to this data is the vector r 
That is 



y - f(x). 



r = [r n ... r n ] 1 , where n = y % - f{x{). 

Our goal is to find / G T such that r is reasonably small. We measure the size of a vector 
by the use of norms, which were explored in Subsection 1.4.1. The most useful norms are the £ p 
norms. For a given p with < p < oo, the £ p norm of r is defined as 

w ' = iS r, J 

For the purposes of approximation, the easiest norm to use is the £ 2 norm: 

Definition 9.1.1 ((Ordinary) Least Squares Best Approximant). The least-squares best 
approximant to a set of data, x, y from a class of functions, T, is the function /* € T that 
minimizes the £ 2 norm of the error. That is, if /* is the least squares best approximant, then 



|y -/*(*) 



mm \\y 



/(*) 
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We will generally assume uniqueness of the minimum. This method is sometimes called the ordinary 
least squares method. It assumes there is no error in the measurement of the data x, and usually 
admits a relatively straightforward solution. 



9.1.2 Linear Least Squares 

We illustrate this definition using the class of linear functions as an example, as this case is reason- 
ably simple. We are assuming 

T = = ax + b\ a,b £R) . 

We can now think of our job as drawing the "best" line through the data. 

By Definition 9.1.1, we want to find the function f(x) = ax + b that minimizes 



\ 1/2 
2 \ 



,i=0 



This is a minimization problem over two variables, x and y. As you may recall from calculus class, 
it suffices to minimize the sum rather than its square root. That is, it suffices to minimize 

n 

^ [axi + b- yi} 2 . 

i=0 

We minimize this function with calculus. We set 



= — = ^2 2x k {ax k + b-y k ) 



a 

k=0 



0= lf = ^2 2 ( aX k + b -Vk) 



k=0 



These are called the normal equations. Believe it or not these are two equations in two unknowns. 
They can be reduced to 

^2 x l a + ^2 x kb = Yx k y k 
^2x k a+ (n + l)b = ^yk 
The solution is found by naive Gaussian Elimination, and is ugly. Let 



'12 



dn 


= £4 


= d 2 i 


= Y, Xk 


d22 


= n+l 


ei 


= £ x k y k 


e2 


= £y fc 



We want to solve 



dna + d 12 b = e\ 
d 2 ia + d 2 2b = e 2 
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Gaussian Elimination produces 

^22^11 — ^12^21 
^22^11 — 6^12^21 

The answer is not so enlightening as the means of finding the solution. 

We should, for a moment, consider whether this is indeed the solution. Our calculations have 
only shown an extrema at this choice of (a, 6); could it not be a maxima? 




9.1.3 Least Squares from Basis Functions 



In many, but not all cases, the class of functions, is the span of a small set of functions. This 
case is simpler to explore and we consider it here. In this case J- can be viewed as a vector space 
over the real numbers. That is, for f,g £ J 7 , and a, (3 G R, then af + (3g £ T , where the function 
af is that function such that (af)(x) = af(x). 

Now let {<7j(^)}™ =0 be a set of m + 1 linearly independent functions, i.e., 



cogo(x) + cigi(x) + . . . + c m g m (x) = Vx 



c = Ci 







Then we say that T is spanned by the functions {gj(x)}™ =Q if 

f( x ) = ^2 c j9j(x) I Cj e R, j = 0, 1, . . . , 



m 



> . 



In this case the functions gj are basis functions for T . Note the basis functions need not be unique: 
a given class of functions will usually have more than one choice of basis functions. 

Example 9.1. The class T = {f(x) = ax + b \ a, b 6 R} is spanned by the two functions go(x) = 1, 
and gi(x) = x. However, it is also spanned by the two functions go(x) = 2x + 1, and g± (x) = x — 1. 

To find the least squares best approximant of T for a given set of data, we minimize the square 
of the £ 2 norm of the error; that is we minimize the function 



0(c o ,ci,...,c m ) = ^2 



k=0 



J2 c j9j{xk) 



Vk 



(9.1) 



Again we set partials to zero and solve 



£2 

fc=0 



Vk 



9i{x k ) 



This can be rearranged to get 

m 

E 

3=0 



J2gj{xk)9i(x k ) 



cj = ^2yk9i{x k ) 



k=0 
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If we now let 

n n 

di 'j = ^9j{x k )gi{x k ), ei = ^y fc ffi(x fe ), 

k=0 k=0 

Then we have reduced the problem to the linear system (again, called the normal equations): 





doi 


d()2 ■ 


dom 




CO 




eo 




du 


di2 ■ 


dim 




Cl 




ei 


d-20 


di\ 


d22 ■ 


d2m 




C2 




e-2 


dmO 


dmX 


d m 2 ■ 


d mm 




Cm 




Cm 



(9.2) 



The choice of the basis functions can affect how easy it is to solve this system. We explore this 
in Section 9.2. Note that we are talking about the basis {gjjjlo, and not exactly about the class 
of functions T . 

For example, consider what would happen if the system of normal equations were diagonal. In 
this case, solving the system would be rather trivial. 



Example Problem 9.2. Consider the case where m = 0, and go{x) = hix. Find the least squares 
approximation of the data 



X 


0.50 


0.75 


1.0 


1.50 


2.0 


2.25 


2.75 


3.0 


y 


-1.187098 


-0.452472 


-0.068077 


0.713938 


1.165234 


1.436975 


1.725919 


1.841422 



Solution: Essentially we are trying to find the c such that clnx best approximates the data. The 
system of equation 9.2 reduces to the 1-D equation: 

[S fc ln 2 x k ] c = Y, k y k lnx k 

For our data, this reduces to: 

4.0960c = 6.9844 

so we find c = 1.7052. The data and the least squares interpolant are shown in Figure 9.1. H 
Example 9.3. Consider the awfully chosen basis functions: 

9o(x) = - l) x 2 - + 1 

gi{x) = x 3 + (|-l) (x 2 + x) +1 

where e is small, around machine precision. 

Suppose the data are given at the nodes xq = 0, x\ = 1,X2 = — 1. We want to set up the normal 
equations, so we compute some of the dij. First we have to evaluate the basis functions at the nodes 
Xj. But this example was rigged to give: 

g (x ) = 1, g {xi) = 0, g {x 2 ) = e 
gi(x ) = 1, gi{x±) = e, gi{x 2 ) = 

After much work we find we want to solve 



1 + e 2 1 




co 




yo + m 


1 1 + e 2 _ 




. c i _ 




_ yo + m _ 
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-2.5 1 1 1 1 1 1 1 1 

0.5 1 1.5 2 2.5 3 

Figure 9.1: The data of Example Problem 9.2 and the least squares interpolant are shown. 

However, the computer would only find this if it had infinite precision. Since it does not, and since 
e is rather small, the computer thinks e 2 = 0, and so tries to solve the system 



" 1 1 " 




CO 






1 1 




c l 




. yo + m _ 



When yi ^ y2, this has no solution. Bummer. 

This kind of thing is common in the method of least squares: the coefficients of the normal 
equations include terms like 

9i(xk)9j(x k ). 

When the gi are small at the nodes x^, these coefficients can get really small, since we are squaring. 

Now we draw a rough sketch of the basis functions. We find they do a pretty poor job of 
discriminating around all the nodes. 

9.2 Orthonormal Bases 

In the previous section we saw that poor choice of basis vectors can lead to numerical problems. 
Roughly speaking, if gi{xk) is small for some i's and /c's, then some dij can have a loss of precision 
when two small quantities are multiplied together and rounded to zero. 

Poor choice of basis vectors can also lead to numerical problems in solution of the normal 
equations, which will be done by Gaussian Elimination. 

Consider the case where T is the class of polynomials of degree no greater than m. For simplicity 
we will assume that all Xi are in the interval [0, 1]. The most obvious choice of basis functions is 
gj(x) = x 3 . This certainly gives a basis for J 7 , but actually a rather poor one. To see why, look 
at the graph of the basis functions in Figure 9.3. The basis functions look too much alike on the 
given interval. 

A better basis for this problem is the set of Chebyshev Polynomials of the first kind, i.e., 
gj{x) = Tj(x), where 



T (x) = 1, Ti(x) = x, T i+ i(x) = 2xTi(x) - Tj_i(x). 
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Figure 9.2: The basis functions 50 0*0 = (f — l) x<1 ~ § x + lj an d 51 (z) = x 3 + (§ — l) (x 2 + x) + 1 
are shown for e = 0.05. Note that around the three nodes 0, 1, —1, these two functions take nearly 
identical values. This can lead to a system of normal equations with no solution. 




Figure 9.3: The polynomials x l for i = 0, 1, ... ,6 are shown on [0, 1]. These polynomials make a 
bad basis because they look so much alike, essentially. 



These polynomials are illustrated in Figure 9.4. 
9.2.1 Alternatives to Normal Equations 

It turns out that the Normal Equations method isn't really so great. We consider other methods. 
First, we define A as the n x m matrix defined by the entries: 



aij = gj(xi), i = 0, 1,... ,n, j = 0,1,... ,m. 
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Figure 9.4: The Chebyshev polynomials Tj(x) for j = 0,1, ... ,6 are shown on [0, 1]. These polyno- 
mials make a better basis for least squares because they are orthogonal under some inner product. 
Basically, they do not look like each other. 



That is 



" 3o Oo) 


31(20) 


32(2o) • 


■ 9m (20) 


30(21) 


91 (21) 


32(2l) • 


• 9m (2l) 


30(22) 


31(22) 


32(22) • 


• 9rn (22) 


30(23) 


31(23) 


32(23) • 


• 9m (23) 


30(24) 


31(^4) 


32(24) • 


• 3m (24) 


. 30 On) 


91 (Xn) 


32 (X n ) ■ 


9m\Xn) 



We write it in this way since we are thinking of the case where n 3> m, so A is "tall." 
After some inspection, we find that the Normal Equations can be written as: 



A T Ac = A T y. 



(9.3) 



Now let the vector c be identified, in the natural way, with a function in J- '. That is c is 
identified with f(x) = YlY=o c j9j( x )- Y° u should now convince yourself that 

Ac = /(*). 

And thus the residual, or error, of this function is r = y — Ac. 

In our least squares theory we attempted to find that c that minimized 

\\y-Ac\\l = {y- Ac) T (y- Ac) 

We can see this as minimizing the Euclidian distance from y to A c. For this reason, we will have 
that the residual is orthogonal to the column space of A, that is we want 



A T r = A T (y- Ac) = 0. 
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This is just the normal equations. We could rewrite this, however, in the following form: find c, r 
such that 



1 A " 




r 




y 


A T 




c 








This is now a system of n + m variables and unknowns, which can be solved by specialized means. 
This is known as the augmented form. We briefly mention that naive Gaussian Elimination is not 
appropriate to solve the augmented form, as it turns out to be equivalent to using the normal 
equations method. 

9.2.2 Ordinary Least Squares in octave/Matlab 

The ordinary least squares best approximant is so fundamental that it is hard-wired into oc- 
tave/Matlab's system solve operator. The general equation we wish to solve is 

Ac = y. 

This system is overdetermined: there are more rows than columns in A, i.e., more equations than 
unknown variables. However, in octave/Matlab, this is solved just like the case where A is square: 

octave : 1> c = A \ y ; 

Thus, unless specifically told otherwise, you should not implement the ordinary least squares solu- 
tion method in octave/Matlab. Under the hood, octave/Matlab is not using the normal equations 
method, and does not suffer the increased condition number of the normal equations method. 

9.3 Orthogonal Least Squares 

The method described in Section 9.1, sometimes referred to as "ordinary least squares," assumes 
that measurement error is found entirely in the dependent variable, and that there is no error in 
the independent variables. 

For example, consider the case of two variables, x and y, which are thought to be related by an 
equation of the form y = mx + b. A number of measurements are made, giving the two sequences 
{xj}™ =1 and {yi}™ = i- Ordinary least squares assumes, that 

yi = mxi + b + €i, 

where ej, the error of the i th measurement, is a random variable. It is usually assumed that 
E[ej] = 0, i.e., that the measurements are "unbiased." Under the further assumption that the 
errors are independent and have the same variance, the ordinary least squares solution is a very 
good one. 1 

However, what if it were the case that 

Hi = m(xi + &) + b + ei, 

i.e., that there is actually error in the measurement of the x{! In this case, the orthogonal least 
squares method is appropriate. 

lr rhis means that the ordinary least squares solution gives an unbiased estimator of m and b, and, moreover, gives 
the estimators with the lowest variance among all linear, a priori estimators. For more details on this, locate the 
Gauss-Markov Theorem in a good statistics textbook. 
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The difference between the two methods is illustrated in Figure 9.5. In Figure 9.5(a), the 
ordinary least squares method is shown; it minimizes the sum of the squared lengths of vertical 
lines from observations to a line, reflecting the assumption of no error in Xj. The orthogonal least 
squares method will minimize the sum of the squared distances from observations to a line, as 
shown in Figure 9.5(b). 




i 1 1 1 1 < 1 1 i 1 1 1 1 < 1 1 

1 2 3 4 5 6 1 2 3 4 5 6 

(a) Ordinary Least Squares (b) Orthogonal Least Squares 



Figure 9.5: The ordinary least squares method, as shown in (a), minimizes the sum of the squared 
lengths of the vertical lines from the observed points to the regression line. The orthogonal least 
squares, shown in (b), minimizes the sum of the squared distances from the points to the line. The 
offsets from points to the regression line in (b) may not look orthogonal due to improper aspect 
ratio of the figure. 

We construct the solution geometrically. Let {x^}^ be the set of observations, expressed as 
vectors in ]R n . The problem is to find the vector n and number d such that the hyperplane n T x = d 
is the plane that minimizes the sum of the squared distances from the points to the plane. We can 
solve this problem using vector calculus. 

The function 

1 m n 2 

J-Y\\n T x^-d . 

is the one to be minimized with respect to n and d. It gives the sum of the squared distances from 
the points to the plane described by n and d. To simplify its computation, we will minimize it 

1 1 1 2 

subject to the constraint that ||n|| 2 = 1. Thus our problem is to find n and d that solve 



min n T £C^ — d 



We can express this as min f(n,d) subject to g(n,d) = 1. The solution of this problem uses the 
Lagrange multipliers technique. 

Let £(n,d,\) = f(n,d) — Xg(n,d). The Lagrange multipliers theory tells us that a necessary 
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condition for n, d to be the solution is that there exists A such that the following hold: 

(§§(n,d,A)=0 
V„£ (n, d, A) = 
K g(n,d) = 1 



Let X be the m x n matrix whose i th row is the vector »W . We can rewrite the Lagrange 
function as 



£ (n, cZ, A) = (Xn — (il m ) (Xn — dl m ) — An T n 

= n T X T Xn - 2dl™ T Xn + d 2 l m T l m - Xn T n 



We solve the necessary condition on 



= — (n, d, A) = 2dl m T l m - 21 m T Xn 



Xn/1 



rn Irr?, 



This essentially tells us that we will zero the first moment of the data around the line. Note that 
this determination of d requires knowledge of n. However, since this is a necessary condition, we 
can plug it into the gradient equation: 



= V n £ (n, d, A) = 2X T Xn - 2d (l m T x) T - 2An = 2 X T Xn - hlL^lH ( 



X 



An 



The middle term is a scalar times a vector, so the multiplication can be commuted, this gives 



X T Xn 



X T 



1 1 1 



-n — Xn 



vTi -i Tv 

X T X - m } m - Al 



1 1 1 



TJ, 



Thus n is an eigenvector of the matrix 



M=X T X- XTl ™ r 1 ™ TX 
1 T 1 

J-m J -m 

The final condition of the three Lagrange conditions is that g{n,d) = n n = 1. Thus if n is a 
minimizer, then it is a unit eigenvector of M. 

It is not clear which eigenvector is the right one, so we might have to check all the eigenvectors. 
However, further work tells us which eigenvector it is. First we rewrite the function to be minimized, 
when the optimal d is used: f(n) = f(n, l TTl T Xn/l TTl T l rn ,). We have 

2 



f(n) = n T X T Xn - 2 



Xn. 



1 T l 

J-m J -m 



Xn + 



lm T Xn 
1 T l 

J-m J -m 



1 T l 



Tv Tv « T X T l m l m T Xn n T X T l m l m T Xnl m T l, 
n X Xn — 2 : — — h - 



n T X T Xn - 2 
n T X T Xn 



1 T l 
n T X T l rT7 l rT7 T Xn n 



+ 



1 T l 1 T l 

T X T l m l m T Xn 



1 T l 



n T X T l m l m T Xn 



n 



T 



x'x 



X T l m l m T X 

I Ti 



n 



n T Mn 
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This sequence of equations only required that the d be the optimal d for the given n, and used the 
fact that a scalar is it's own transpose, thus l m T Xn = ra T X T l m . 

Now if n is a unit eigenvector of M, with corresponding eigenvalue A, then n T Mn = A. Note 
that because f(n,d) is defined as w T w, for some vector w, then the matrix M must be positive 
definite, i.e., A must be positive. However, we want to minimize A. Thus we take n to be the unit 
eigenvector associated with the smallest eigenvalue. 

Note that, by roundoff, you may compute that M has some negative eigenvalues, though they 
are small in magnitude. Thus one should select, as n, the unit eigenvector associated with the 
smallest eigenvalue in absolute value. 

Example Problem 9.4. Find the equation of the line which best approximates the 2D data 
{(xi,yi)}™ =1 , by the orthogonal least squares method. 
Solution: In this case the matrix X is 



Xl 


Vi 


X2 


2/2 


^3 


2/3 



Thus we have that 



m = x t x _x^o^x 



1 

m 



— m 



x xy 
xy y 2 



E x i E Xiy, 

E x iVt E vl 

E x i E x iV- 
E x iVi E yf 

E x\ — rax 1 E x iy-i — mxy 
E Xiyi - mxy E yf ~ m y 2 
2a (3 
13 2 7 



(E Xi) E Xi E yi 



E Xi E yi (E yiY 



(9.4) 



where we use x, y, to mean, respectively, the mean of the x- and y-values. The characteristic 
polynomial of the matrix, whose roots are the eigenvalues of the matrix is 



p(A) = det 



2a- A 
f3 2 7 -A 



= 4a 7 -/3 2 -2(a + 7 )A + A 2 . 



The roots of this polynomial are given by the quadratic equation 
2(q + 7) ± ^4(a + 7 ) 2 -4(4a 7 -/3 2 ) 



A± 



(a + 7)± V(«-7) +/3 2 



We want the smaller eigenvalue, A_ = (a + 7) — ^/ (a — r y) 2 + j3 2 . The associated eigenvector is in 
the null space of the matrix 



= 



2a - A_ 
P 2 7 - A_ 



v = 



2a - A_ /3 
P 2 7 -A. 



-A: 
1 



/?- k(2a- A_) 
-A;/? + 27- A_ 
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This is solved by 

2 7 - A_ _ (7 " ") + \Jh~af + P 2 



k = 



(9.5) 



YsiVi -xj)-m (y 2 - x 2 ) + y {Y,{yf ~ tf) -m(y 2 - x 2 )) 2 + 4 (]T Xi yi - mxy) 2 

2E Xiyi ~ mxy 
Thus the best plane through the data is of the form 

— kx + ly = d or y = kx + d 

Note that the eigenvector, v, that we used is not a unit vector. However it is parallel to the unit 
vector we want, and can still use the regular formula to compute the optimal d. 



d — l m T Xn/l T7T T l T7T — — lm. T X 

m 



-k 
1 



[ x y 



-k 
1 



y — kx 



Thus the optimal line through the data is 

y = kx + d= k(x — x)+y or y — y = k(x — x) , 
with k defined in equation 9.5. 



9.3.1 Computing the Orthogonal Least Squares Approximant 

Just as the Normal Equations equation 9.3 define the solution to the ordinary least squares ap- 
proximant, but are not normally used to solve the problem, so too is the above described means of 
computing the orthogonal least squares not a good idea. 
First we define 



W = X 



1 7n 1 rrv, X 
I Ti • 



Now note that 



W W - i X - lml ™ TX ) (x - lml r TX ) X X - 2 



1 1 1 



1 1 1 



X T lmlm T X (lmlm^X) l^l^^X 
lm T lrTi (l m T lm) 2 

V T V r) X T l T 7alTTi T X X T l Tn ,(l rn T l rn ,)l TT1 T X j X T l m l m T X X T l Tn ,l T17 , T X 
— A A — Z =p 1 — ^ — A A — z =p 1 =p 



1 T l 



(1 T i V 



1 T l 
J-tti - L m. 



1 T l 



Y T 1 1 Tv 

= X T X- ™ m = M. 
1 T 1 



We now need some linear algebra magic: 
Definition 9.3.1. A square matrix U is called unitary if its inverse is its transpose, i.e., 

U T U = I = UU T . 

From the first equation, we see that the columns of U form a collection of orthonormal vectors. 
The second equation tells that the rows of U are also such a collection. 
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Definition 9.3.2. Every m x n matrix B can be decomposed as 

B = UZV T , 

where U and V are unitary matrices, U is m x m, V is n x n, and Y. is m x n, and has nonzero 
elements only on its diagonal. The values of the diagonal of T. are the singular values of B. The 
column vectors of V span the row space of B, and the column vectors of U contain the column space 
of B. 

The singular value decomposition generalizes the notion of eigenvalues and eigenvalues to the 
case of nonsquare matrices. Let be the i th column vector of U, the i th column of V, and 
let <7j be the i th element of the diagonal of Z. Then if x = J2i otiv( l \ we have 



Bx = UZV 



• CXr. 



G\OL\ 
O20L2 

























OiOLiU 



(0 



Thus t>W and act as an "eigenvector pair" with "eigenvalue" cij. 

The singular value decomposition is computed by the octave/Matlab command [U,S,V] = 
svd(B). The decomposition is computed such that the diagonal elements of S, the singular values, 
are positive, and decreasing with increasing index. 

Now we can use the singular value decomposition on W: 



W = UZV 



T 



Thus 

M = W T W = VZ T U T UZV T = VZ T IZV T = VS 2 V T , 

where S 2 is the n x n diagonal matrix whose i th diagonal is of. Thus we see that the solution to 
the orthogonal least squares problem, the eigenvector associated with the smallest eigenvalue of M, 
is the column vector of V associated with the smallest, in absolute value, singular value of W. In 
octave/Matlab, this is V(: ,n). 

This may not appear to be a computational savings, but if M is computed up front, and its 
eigenvectors are computed, there is loss of precision in its smallest eigenvalues which might lead us 
to choose the wrong normal direction (especially if there is more than one!). On the other hand, 
M is a n x n matrix, whereas W is m x n. If storage is at a premium, and the method is being 
reapplied with new observations, it may be preferrable to keep M, rather than keeping W. That is, 
it may be necessary to trade conditioning for space. 



9.3.2 Principal Component Analysis 

The above analysis leads us to the topic of Principal Component Analysis. Suppose the m x n 
matrix X represents m noisy observations of some process, where each observation is a vector in 
W 1 . Suppose the observations nearly lie in some /c-dimensional subspace of R n , for k < n. How can 
you find an approximate value of k, and how do you find the /c-dimensional subspace? 
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As in the previous subsection, let 

W = X - hn^L^ = X - l m X. where X = l m T X/l m T l m . 

Now note that X is the mean of the m observations, as a row vector in M n . Under the assumption 
that the noise is approximately the same for each observation, we should think that X is likely 
to be in or very close to the fc-dimensional subspace. Then each row of W should be like a vector 
which is contained in the subspace. If we take some vector v and multiply Wd, we expect the 
resultant vector to have small norm if v is not in the /c-dimensional subspace, and to have a large 
norm if it is in the subspace. 

You should now convince yourself that this is related to the singular value decomposition of W. 
Decomposing v = Y2i a i v ^\ we have \Nv = dionu^, and thus product vector has small norm if 
the a.% associated with large cij are small, and large norm otherwise. That is the principal directions 
of the "best" fc-dimensional subspace associated with X are the k columns of V associated with the 
largest singular values of W. 




0.3 1 1 1 1 1 1 1 1 1 

5 10 15 20 25 30 35 40 

Figure 9.6: Noisy 5-dimensional data lurking in M are treated to Principal Component Analysis. 
This figure shows the normalized cumulative sum of squared singular values of W. There is a sharp 
"elbow" at k = 5, which indicates the data is actually 5-dimensional. 

If k is unknown a priori, a good value can be obtained by eyeing a graph of the cumulative 
sum of squared singular values of W, and selecting the k that makes this appropriately large. This 
is illustrated in Figure 9.6, with data generated by the following octave/Matlab code: 

octave :1> X = rand(150,5) ; 

octave:2> sawX = (X * randn(5,40)) + kron(ones(150, 1) ,randn(l,40)) ; 
octave :3> sawX += 0.1 * randn(150,40) ; 
octave:4> wun = ones(150,l); 
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octave :5> W = sawX - (wun * wun' * sawX) / (wun' * wun) ; 

octave:6> [U,S,V] = svd(W) ; 

octave :7> eigs = diag(S'*S); 

octave:8> cseig = cumsum(eigs) ./ sum(eigs) ; 

octave :9> plot ( (1 : 40) ' , cseig, 'xr-@ ; cum. sum sqrd. sing, vals.;'); 
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Exercises 

(9.1) How do you know that the choice of constants Cj in our least squares analysis actually find a 
minimum of equation 9.1, and not, say, a maximum? 

(9.2) Our "linear" least squares might be better called the "affine" least squares. In this exercise 
you will find the best linear function which approximates a set of data. That is, find the 
function f(x) = cx which is the ordinary least squares best approximant to the given data 



X 


x 


Xi 






y 


yo 


Vi 




Vn 



(9.3) Find the constant that best approximates, in the ordinary least squares sense, the given data 
x, y. (Hint: you can use equation 9.2 or equation 9.3 using a single basis function go(x) = 1.) 
Do the X{ values affect your answer? 

(9.4) Find the function f(x) = c which best approximates, in the ordinary least squares sense, the 
data 



X 


l 


-2 


5 


y 


l 


-2 


4 



(9.5) Find the function ax + b that best approximates the data 



X 





-l 


2 


y 





l 


-l 



Use the ordinary least squares method. 

(9.6) Find the function ax + b that best approximates the data of the previous problem using the 
orthogonal least squares method. 

(9.7) Find the function ax 2 + b that best approximates the data 



X 





1 


2 


y 


0.3 


0.1 


0.5 



(9.8) Prove that the least squares solution is the solution of the normal equations, equation 9.3, 
and thus takes the form 

c= ^A T A^ 1 A T y. 

(9.9) For ordinary least squares regression to the line y = mx + b, express the approximate m 
in terms of the a, j3, and 7 parameters from equation 9.4. Compare this to that found in 
equation 9.5. 

(9.10) Any symmetric positive definite matrix, G can be used to define a norm in the following 
way: 

\\v\\ G =df V v T Gv 

The weighted least squares method for G and data A and y, finds the c that solves 

min || Ac — y||^ 

Prove that the normal equations form of the solution of this problem is 

A T GAc = A T Gy. 
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(9.11) (Continuation) Let the symmetric positive definite matrix G have Cholesky Factorization 
G = LL T , where L is a lower triangular matrix. Show that the solution to the weighted least 
squares method is the c that is the ordinary (unweighted) least squares best approximate 
solution to the problem 



L T Ac = L T i 



(9.12) The r 2 statistic is often used to describe the "goodness-of-fit" of the ordinary least squares 
approximation to data. It is defined as 



r 



2 l|Ac|| 2 

II l|2 ' 
Hi/112 



where c is the approximant A T A 1 A T y. 
(a) Prove that we also have 



r 



1 _ 1 _ l|Ac — y||| 

II l|2 



(b) Prove that the r 2 parameter is "scale- invariant," that is, if we change the units of our 
data, the parameter is unchanged (although the value of c may change.) 

(c) Devise a similar statistic for the orthogonal least squares approximation which is scale- 
invariant and rotation-invariant. 

(9.13) As proved in Example Problem 9.4, the orthogonal least squares best approximate line for 
data {{xi,yi)}™ =1 goes through the point (x,y). Does the same thing hold for the ordinary 
least squares best approximate line? 

(9.14) Find the constant c such that f(x) = In (cx) best approximates, in the least squares sense, 
the given data 



X 


xo 


X\ 




Xn 


y 


yo 


yi 




Vn 



(Hint: You cannot use basis functions and equation 9.2 to solve this. You must use the 
Definition 9.1.1.) The geometric mean of the numbers a\, 0,2, . . . ,a n is defined as (Y\ a-i) l ^ n ■ 
How does your answer relate to the geometric mean? 

(9.15) Write code to find the function ae x + be~ x that best interpolates, in the least squares sense, 
a set of data {(xi, yi)Y~^. Your m-file should have header line like: 

function [a,b] = lsqrexp(xys) 

where xys is the (n + 1) x 2 matrix whose rows are the n + 1 tuples (xi,yi). Use the normal 
equations method. This should reduce the problem to solving a 2 x 2 linear system; let 
octave/Matlab solve that linear system for you. (Hint: To find x such that Mx = b, use x = 
M \ b.) 

(a) What do you get when you try the following? 
octave :1> xs = [-1 -0.5 0.5 1]'; 

octave : 2> ys = [1.194 0.430 0.103 0.322 1.034]'; 
octave :3> xys = [xs ys] ; 
octave :4> [a,b] = lsqrexp(xys) 

(b) Try the following: 

octave :5> xys = xys = [-0.00001 0.3;0 0.6; 0.00001 0.7]; 
octave :6> [a,b] = lsqrexp(xys) 

Does something unexpected happen? Why? 

(9.16) Write code to find the plane n T x^ = d that best interpolates, in the least squares sense, a 
set of data {x^Y i= ™- Your m-file should have header line like: 
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function [n,d] = orthlsq(X) 

where X is the m x n matrix whose rows are the m tuples (xi, X2, ■ ■ • , x n ). Use the oc- 
tave/Matlab eig function, which returns the eigenvalues and vectors of a matrix. 

(a) Create artificially planar data by the following: 
octave:l> m = 100;n = 5;epsi = 0.1;offs = 2; 
octave :2> Xsub = rand(m, (n-1) ) ; 

octave :3> T = rand (n-1 ,n) ; 

octave :4> X = Xsub * T; 

octave :5> X += offs .* ones (size (X) ) ; 

octave :6> X += epsi .* randn(size(X) ) ; 

Now X contains data which are approximately on a (n — l)-dimensional hyperplane in 
W l . The data, when centered, lay in the rowspace of T. Now test your code. 
octave:7> [n,d] = orthlsq(X) ; 

You should have that the vector n is orthogonal to the rows of T. Check this: 
octave :8> disp(norm(T * n)); 

What do you get? Try again with epsi very small (1 x 10 ~ 5 ) or zero. 

(b) Compare the orthogonal least squares technique to the ordinary least squares in a two 
dimensional setting. Consider x, y data nearly on the line y = mx + b: 

octave :1> obs = 100;xeps = 1.0;yeps = 1.0 ;m = 0.375;b = 1; 
octave :2> Xtrue = randn(obs , 1) ; Ytrue = m .* Xtrue .+ b; 
octave :3> XSaw = Xtrue + xeps * randn(size (Xtrue) ) ; 
octave :4> YSaw = Ytrue + yeps * randn (size (Ytrue) ) ; 
octave :5> lsmb = [XSaw, ones (obs, 1)] \ YSaw; 
octave :6> Ism = lsmb(l); lsb = lsmb (2) ; 
octave :7> [n,d] = orthlsq( [XSaw, YSaw]); 
octave :8> orthm = -n(l) / n(2);orthb = d / n(2) ; 

Compare the estimated coefficients from both solutions. Plot the two regression lines, 
octave :9> hold off; 

octave: 10> plot (XSaw, YSaw, "*; observed data;"); 
octave :11> hold on; 

octave: 12> plot (XSaw, m*XSaw .+ b,"r-;true line;"); 

octave :13> plot (XSaw, lsm*XSaw .+ lsb, "g-; ordinary least squares;"); 
octave :14> plot (XSaw, orthm*XSaw .+ orthb, "b-; orthogonal least squares;"); 
octave :15> hold off; 

Try this again with different values of xeps and yeps, which control the sizes of the 
errors in the x and y dimensions. Try with xeps = 0.1, yeps = 1.0, and with with 
xeps = 1.0, yeps = 0.1. Under which situations does orthogonal least squares outperform 
ordinary least squares? Do they give equally erroneous results in some situations? Can 
you think of a way of "reweighting" channel data to make the orthogonal least squares 
method more robust in cases of unequal error variance? 



Chapter 10 

Ordinary Differential Equations 



Ordinary Differential Equations, or ODEs, are used in approximating some physical systems. Some 
classical uses include simulation of the growth of populations and trajectory of a particle. They are 
usually easier to solve or approximate than the more general Partial Differential Equations, PDEs, 
which can contain partial derivatives. 

Much of the background material of this chapter should be familiar to the student of calculus. 
We will focus more on the approximation of solutions, than on analytical solutions. For more 
background on ODEs, see any of the standard texts, e.g., [8]. 



10.1 Elementary Methods 

A one-dimensional ODE is posed as follows: find some function x(t) such that 



In calculus classes, you probably studied the case where f(t,x(t)) is independent of its second 
variable. For example the ODE given by 



A more general case is when f(t,x) is separable, that is f(t,x) = g(t)h{x) for some functions 
g, h. You may recall that the solution to such a separable ODE usually involved the following 
equation: 



Finding an analytic solution can be considerably more difficult in the general case, thus we turn 
to approximation schemes. 




(10.1) 




has solution x(t) = — ^t 2 + K, where K is chosen such that x(a) = c. 
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10.1.1 Integration and 'Stepping' 

We attempt to solve the ODE by integrating both sides. That is 

*^± = f(t,x(t)), yields 

rt+h r t+h 



/t+n rt+n 
dx = f(r,x(r))dr, thus 

/t+h 
f{r,x{r))dr. (10.2) 



If we can approximate the integral then we have a way of 'stepping' from t to t + h, i.e., if we have 
a good approximate of x(t) we can approximate x(t + h). Note that this means that all the work 
we have put into approximating integrals can be used to approximate the solution of ODEs. 

Using stepping schemes we can approximate x(tfi na i) given x(ti n i t i a i) by taking a number of 
steps. For simplicity we usually use the same step size, h, for each step, though nothing precludes 
us from varying the step size. 

If we apply the left-hand rectangle rule for approximating integrals to equation 10.2, we get 



x(t + h)^x(t) + hf(t,x(t)). (10.3) 



This is called Euler's Method. Note this is essentially the same as the "forward difference" approx- 
imation we made to the derivative, equation 7.1. 
Trapezoid rule gives 

x(t + /i) » x{t) + ^ [f(t, x{t)) + f(t + h, x(t + h))\ . 

But we cannot evaluate this exactly, since x(t + h) appears on both sides of the equation. And it is 
embedded on the right hand side. Bummer. But if we could approximate it, say by using Euler's 
Method, maybe the formula would work. This is the idea behind the Runge-Kutta Method (see 
Section 10.2). 



10.1.2 Taylor's Series Methods 

We see that our more accurate integral approximations will be useless since they require information 
we do not know, i.e., evaluations of f(t,x) for yet unknown x values. Thus we fall back on 
Taylor's Theorem (Theorem 1.1). We can also view this as using integral approximations where all 
information comes from the left-hand endpoint. 

By Taylor's theorem, if x has at least m + 1 continuous derivatives on the interval in question, 
we can write 

x (t + h)= x{t) + hx'{t) + \h 2 x"{t) + ... + —h m x^ m \t) + - — l —-h m+1 x^ m+1 \T), 

where r is between t and t + h. Since r is essentially unknown, our best calculable approximation 
to this expression is 

x{t + h)r* x(t) + hx'{t) + \h 2 x"{t) + ... + —h m x^ (t). 

This approximate solution to the ODE is called a Taylor's series method of order m. We also say 
that this approximation is a truncation of the Taylor's series. The difference between the actual 
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value of x(t + h), and this approximation is called the truncation error. The truncation error exists 
independently from any error in computing the approximation, called roundoff error. We discuss 
this more in Subsection 10.1.5. 



10.1.3 Euler's Method 

When m = 1 we recover Euler's Method: 



x{t + h) = x(t) + hx'(t) = x{t) + hf(t, x{t)). 



Example 10.1. Consider the ODE 



dt ~ x ' 

z(0) = 1. 



The actual solution is x{t) = e*. Euler's Method will underestimate x(t) because the curvature of 
the actual x(t) is positive, and thus the function is always above its linearization. In Figure 10.1, 
we see that eventually the Euler's Method approximation is very poor. 



fculer approxir^atioij 




1 1 1 1 1 1 1 

0.5 1 1.5 2 2.5 3 

Figure 10.1: Euler's Method applied to approximate x' = x, x(0) = 1. Each step proceeds on 
the linearization of the curve, and is thus an underestimate, as the actual solution, x(t) = e has 
positive curvature. Here step size h = 0.3 is used. 



10.1.4 Higher Order Methods 

Higher order methods are derived by taking more terms of the Taylor's series expansion. Note 
however, that they will require information about the higher derivatives of x(t), and thus may not 
be appropriate for the case where f(t,x) is given as a "black box" function. 
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Example Problem 10.2. Derive the Taylor's series method for m = 3 for the ODE 

dx( fil — T I p x 
At ' 

x(0) = 1. 

Solution: By simple calculus: 

x'(t) = x + e x 

x "(t) = x' + e x x' 

x"'(t) = x"(l + e x ) + x'e x x' 

We can then write our step like a program: 

x'(t) <- x(t) + e x(t) 

z"(t) «-x'(t) (l + e x ^ 

x »>(t) <_ x "( t ) (l + e -W) + V« 

x(/ + /i) « x(t) + hx'(t) + J/iV(t) + -/i 3 x'"(t) 

— ■ D 



10.1.5 Errors 



There are two types of errors: truncation and roundoff. Truncation error is the error of truncating 
the Taylor's series after the m th term. You can see that the truncation error is O (/i m+1 ). Roundoff 
error occurs because of the limited precision of computer arithmetic. 

While it is possible these two kinds of error could cancel each other, given our pessimistic 
worldview, we usually assume they are additive. We also expect some tradeoff between these two 
errors as h varies. As h is made smaller, the truncation error presumably decreases, while the 
roundoff error increases. This results in an optimal /i-value for a given method. See Figure 10.2. 
Unfortunately, the optimal /i-value will usually depend on the given ODE and the approximation 
method. For an ODE with unknown solution, it is generally impossible to predict the optimal h. We 
have already observed this kind of behaviour, when analyzing approximate derivatives-c/. Example 
Problem 7.5. 

Both of these kinds of error can accumulate: Since we step from x(t) to x(t + h), an error in 
the value of x(t) can cause a greater error in the value of x(t + h). It turns out that for some ODEs 
and some methods, this does not happen. This leads to the topic of stability. 

10.1.6 Stability 

Suppose you had an exact method of solving an ODE, but had the wrong data. That is, you were 
trying to solve 

x(a) = c, 

but instead fed the wrong data into the computer, which then solved exactly the ODE 

^ =/&*(*))> 
x(a) = c + e. 
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truncation _ _ _ _ 
roundoff 




1e-07 1e-06 1e-05 0.0001 0.001 0.01 

Figure 10.2: The truncation, roundoff, and total error for a fictional ODE and approximation 
method are shown. The total error is apparently minimized for an h value of approximately 0.0002. 
Note that the truncation error appears to be O (h), thus the approximation could be Euler's Method. 
In reality, we expect the roundoff error to be more "bumpy," as seen in Figure 7.1. 



This solution can diverge from the correct one because of the different starting conditions. We 
illustrate this with the usual example. 

Example 10.3. Consider the ODE 

dx(t) 

— ; — = x - 
dt 

This has solution x(t) = x(0)e . The curves for different starting conditions diverge, as in Fig- 
ure 10.3(a). 

However, if we instead consider the ODE 

dx(t) 

which has solution x(t) = x(0)e , we see that differences in the initial conditions become immate- 
rial, as in Figure 10.3(b). 

Thus the latter ODE exhibits stability: roundoff and truncation errors accrued at a given step 
will become irrelevant as more steps are taken. The former ODE exhibits the opposite behavior- 
accrued errors will be amplified. 
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Figure 10.3: Exponential functions with different starting conditions: in (a), the functions (l + e)e 
are shown, while in (b), (l+e)e~* are shown. The latter exhibits stability, while the former exhibits 
instability. 
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It turns out there is a simple test for stability. If f x > 6 for some positive 8, then the ODE 
is instable; if f x < —5 for a positive 5, the equation is stable. Some ODEs fall through this test, 
however. 

We can prove this without too much pain. Define x(t, s) to be the solution to 

^ = f(t,x(t)), 
x(a) = s, 



for t > a. 

Then instability means that 

d 



lim 

t^oo 



oo. 



Think about this in terms of the curves for our simple example x' = x. 
If we take the derivative of the ODE we get 

d 

— x(t,s) = f(t,x(t,s)), 
^| x(M) = ^ /(t ' x(M)) ' Z ' e " 

^"^U(t, S) = f t (t, S(t))^ + f x (t, X(t, S)) 9X{ 



os at as as 

By the independence of t, s the first part of the RHS is zero, so 

||x(M) = / ,(t,x(M))^M. 

We use continuity of x to switch the orders of differentiation: 

9 9, \ , , , u 9 , . 
— — x{t,s) = f x (t,x{t,s))—x{t,s). 

Defining u(t) = -^x(t,s), and q(t) = f x (t,x(t,s)), we have 

^u(t)=q(t)u(t). 

The solution is u(t) = ce^\ where 

Q(t)= f q{r)dr. 

J a 

If q(r) = f x (r,x(r,s)) > 5 > 0, then lim Q(t) = oo, and so limu(i) = oo. But note that u(t) = 
■§^x(t, s), and thus we have instability. 

Similarly if q(r) < — 5 < 0, then limQ(t) = — oo, and so limu(t) = 0, giving stability. 

10.1.7 Backwards Euler's Method 

Euler's Method is the most basic numerical technique for solving ODEs. However, it may suffer 
from instability, which may make the method inappropriate or impractical. However, it can be 
recast into a method which may have superior stability at the cost of limited applicability. 
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This method, the so-called Backwards Euler's Method, is a stepping method: given a reasonable 
approximation to x(t), it calculates an approximation to x(t + h). As usual, Taylor's Theorem is 
the starting point. Letting tf = t + h, Taylor's Theorem states that 

x(t f -h) = x(t f ) + (-h)x'(t f ) + ^-x"(t f ) + ..., 

for an analytic function. Assuming that x has two continuous derivatives on the interval in question, 
the second order term is truncated to give 

x(t) =x{t + h-h)= x(t f - /i) ~ x(t f ) + (-h)x'(t } ) 
and so, x(t + h) « x(t) + hx'(t + h) 



The complication with applying this method is that x'(t+h) may not be computable if x(t+h) is 
not known, thus cannot be used to calculate an approximation for x(t+h). Since this approximation 
may contain x(t + h) on both sides, we call this method an implicit method. In contrast, Euler's 
Method, and the Runge-Kutta methods in the following section are explicit methods: they can be 
used to directly compute an approximate step. We make this clear by examples. 



Example 10.4. Consider the ODE 



dx(t) 

~St = cosx > 
x(0) = 2vr. 



In this case, if we wish to approximate x(0 + h) using Backwards Euler's Method, we have to find 
x{h) such that x(h) = x(0) + cosx(h). This is equivalent to finding a root of the equation 

g(y) = 2vr + cosy - y = 

The function g(y) is nonincreasing (g'(y) is nonpositive), and has a zero between 6 and 8, as seen 
in Figure 10.4. The techniques of Chapter 4 are needed to find the zero of this function. 




Figure 10.4: The nonincreasing function g(y) = 2tt + cosy — y is shown. 
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Example 10.5. Consider the ODE from Example 10.1: 

dx(t) 



dt 

x(0) 



■■ X, 
1. 



The actual solution is x(t) = e t . Euler's Method underestimates x(t), as shown in Figure 10.1. 
Using Backwards Euler's Method, gives the approximation: 

x(t + h) = x(t) + hx(t + h), 
x(t) 



x(t + h) 



1 - h 



Using Backwards Euler's Method gives an overestimate, as shown in Figure 10.5. This is because 
each step proceeds on the tangent line to a curve ke t to a point on that curve. Generally Backwards 
Euler's Method is preferred over vanilla Euler's Method because it gives equivalent stability for 
larger stepsize h. This effect is not evident in this case. 




Figure 10.5: Backwards Euler's Method applied to approximate x' = x, x(0) = 1. The approxi- 
mation is an overestimate, as each step is to a point on a curve ke t , through the tangent at that 
point. Compare this figure to Figure 10.1, which shows the same problem approximated by Euler's 
Method, with the same stepsize, h = 0.3. 



10.2 Runge-Kutta Methods 

Recall the ODE problem: find some x(t) such that 

f(t,x(t)), 



dx(t) 
dt 
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where /, a, c are given. 

Recall that the Taylor's series method has problems: either you are using the first order method 
(Euler's Method), which suffers inaccuracy, or you are using a higher order method which requires 
evaluations of higher order derivatives of x(t), i.e., x"(t), x'"(t), etc. This makes these methods less 
useful for the general setting of / being a blackbox function. The Runge-Kutta Methods (don't 
ask me how it's pronounced) seek to resolve this. 

10.2.1 Taylor's Series Redux 

We fall back on Taylor's series, in this case the 2-dimensional version: 

°° 1 / 8 8\ i 
f(* + h,y + k) = ^-(h- + k-) f(x,y). 

The thing in the middle is an operator on f(x,y). The partial derivative operators are interpreted 
exactly as if they were algebraic terms, that is: 

K'Tx +k Ty) I(X ' V) 

There is a truncated version of Taylors series: 

/(, + h,y + k) - g g [k Tx + k-j /(«,,) 4- (» s + k-j /(«,«. 

where (x, y) is a point on the line segment with endpoints (x, y) and (x + h, y + k). 
For practice, we will write this for n = 1: 

f(x + h,y + k) = f(x, y) + hf x (x, y) + kf y (x, y) + - (h 2 f xx (x, y) + 2hkf xy (x, y) + k 2 f yy (x, y)) 

10.2.2 Deriving the Runge-Kutta Methods 

We now introduce the Runge-Kutta Method for stepping from x(t) to x(t + h). We suppose that 
a, f3,w\,W2 are fixed constants. We then compute: 

K X i-hf(t,x) 

K 2 <- hf(t + ah,x + pKi}. 

Then we approximate: 

x(i + h)<r- x(t) + wi^Ti + w 2 K 2 . 

We now examine the "proper" choices of the constants. Note that w\ = l,w 2 = corresponds 
to Euler's Method, and does not require computation of K 2 . We will pick another choice. 



= f(x,y), 



h df(x,y) h df(x,y) 



dx 



dy 



h 2d ^m^ + 2 hk ^im + k 2d2f ^ 



dx 2 



dxdy 



dy 2 
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Also notice that the definition of K 2 should be related to Taylor's theorem in two dimensions. 
Let's look at it: 

K 2 /h = f(t + ah,x + PK 1 ) 

= f(t + ah,x + (3hf(t,x)) 

= f + ahf t + /3hff x + \ (a 2 h 2 f tt (i, x) + af3fh 2 f tx (i, x) + 2 h 2 f 2 f x x{i, x)) . 
Now reconsider our step: 

x(t + h)= x(t) + Wl K x + w 2 K 2 

= x(t) + Wl hf + w 2 hf + w 2 ah 2 f t + w 2 (3h 2 ff x + O (h 3 ) . 
= x(t) + ( Wl + w 2 ) hx'(t) + w 2 h 2 (aft + pff x ) + O (h 3 ) . 
= x(t) + ( Wl + w 2 ) hx'(t) + w 2 h 2 { a f t + Pff x ) + O (h 3 ) . 
If we happen to choose our constants such that 



wi + w 2 = 1, aw 2 = - = 0w 2 , 



then we get 



x(t + h) = x(t) + hx'(t) + -h 2 {f t + ff x ) + O (h 3 ) 

= x(t) + hx'(t) + hi 2 x"{t) + O (h 3 ) , 

i.e., our choice of the constants makes the approximate x(t + h) good up to a O (/i 3 ) term, because 
we end up with the Taylor's series expansion up to that term. cool. 

The usual choice of constants is a = (5 = 1, w\ = w 2 = \. This gives the second order Runge- 
Kutta Method : 

x(t + h)^ x(t) + -f{t, x) + -/(/ + h, x + hf(t, x)). 
This can be written (and evaluated) as 



K!<-hf(t,x) 

K 2 *- hf(t + h,x + Ki) 

x(t + h)<- x(t) + - (Kx + K 2 ) 



Another choice is a = (3 = 2/3, w\ = l/4,w 2 = 3/4. This gives 



h 3/i 
x (t + h)^x(t) + -f(t,x) + —f[l 



(10.4) 



The Runge-Kutta Method of order two has error term O {h 3 ). Sometimes this is not enough 
and higher-order Runge-Kutta Methods are used. The next Runge-Kutta Method is the order four 
method: 





<-hf(t,x) 


K 2 


«- hf{t + h/2,x + K 1 /2) 


K 3 


^hf(t + h/2,x + K 2 /2) 


K 4 


^hf(t + h,x + K 3 ) 


x(t + h) 


x(t) + \(Ki + 2K 2 + 2K 3 + K 4 ) . 



(10.5) 
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This method has order O (/i 5 ) . (See http : //mathworld. wolf ram. com/Runge-KuttaMethod.html) 
The Runge-Kutta Method can be extrapolated to even higher orders. However, the number 

of function evaluations grows faster than the accuracy of the method. Thus the methods of order 

higher than four are normally not used. 

We already saw that w\ = 1, W2 = corresponds to Euler's Method. We consider now the case 

wi = 0, W2 = 1. The method becomes 

X (t + h )^ x ( t ) + hf(t+^,x+ lf(t, s)) . 

This is called the modified Euler's Method . Note this gives a different value than if Euler's 
Method was applied twice with step size h/2. 

10.2.3 Examples 

Example 10.6. Consider the ODE: 

f x' = (te) s -(f) 2 
\ x(l) = 1 

Use h = 0.1 to compute x(l.l) using both Taylor's Series Methods and Runge-Kutta methods 
of order 2. 



10.3 Systems of ODEs 

Recall the regular ODE problem: find some x(t) such that 

^ = f(t,x(t)), 
x(a) = c, 

where /, a, c are given. 

Sometimes the physical systems we are considering are more complex. For example, we might 
be interested in the system of ODEs: 

^l = f(t,x(t),y(t)), 

*M=g(t,x(t),y(t)), 
x(a) = c, 
y(a) = d. 

Example 10.7. Consider the following system of ODEs: 



f dx(t) 
dt 



t 2 - X 



dt — y t y l, 
x(0) = 1, 

{ y(o) = o. 

You should immediately notice that this is not a system at all, but rather a collection of two ODEs: 



dx(t) 
dt 

x(0) = 1, 



t 2 - X 
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and 



Ml -v 2 + v 
dt — y ^ y 

2/(0) = 0. 



These two ODEs can be solved separately. We call such a system uncoupled. In an uncoupled 
system the function f(t,x,y) is independent of y, and g(t,x,y) is independent of x. A system 
which is not uncoupled, is, of course, coupled. We will not consider uncoupled systems. 



10.3.1 Larger Systems 

There is no need to stop at two functions. We may imagine we have to solve the following problem: 
find x\ (t) , X2 (t) , . . . , x n {t) such that 

f ^ = fi(t, Xl (t),x 2 (t),...,x n (t)), 
<^ = f 2 (t,x 1 (t),x 2 (t),...,x n (t)), 

^■ = fn{t,X 1 {t),X 2 (t),...,X n (t)), 

xi(a) = ci, x 2 (a) = c 2 , . . .,x n (a) = c n . 

The idea is to not be afraid of the notation, write everything as vectors, then do exactly the 
same thing as for the one dimensional case! That's right, there is nothing new but notation: 
Let 





Xl 




- x' } - 




' fi ' 




Cl 


X = 


X2 


, x' = 


x 2 


, F = 


h 


, c = 


C2 








. X n . 




fn 




Cn 



Then we want to find X such that 



X'(t) = F(t,X(t)) 
X (a) = C. 



(10.6) 



Compare this to equation 10.1. 

10.3.2 Recasting Single ODE Methods 

We will consider stepping methods for solving higher order ODEs. That is, we know X(t), and 
want to find X (t + h) . Most of the methods we looked at for solving one dimensional ODEs can 
be rewritten to solve higher dimensional problems. 

For example, Euler's Method, equation 10.3 can be written as 



X(t + h)^ X(t) + hF (t, X(t)) . 



As in ID, this method is derived from approximating X by its linearization. 
In fact, our general strategy of using Taylor's Theorem also carries through without change. 
That is we can use the k th order method to step as follows: 

X (t + h)^ x(t) + hx'(t) + |x"(t) + . . . + ^x {h) (t). 
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Additionally we can write the Runge-Kutta Methods as follows: 



Order two: 





<- hF (t,X) 


K-2 


^hF(t + h,X + K 1 ) 


X(t + h) 


^X(t) + ±(K 1 + K 2 ). 



Order four: 



K x <- hF(t,X) 
K 2 <— hF (^ t + h,X + ^K 1 

K s <— hF ( t +±h,X + ^K 2 
K 4 ^hF(t + h,X + K 3 ) 
X(t + h)^ x(t) + i (JTi + 2K 2 + 2K 3 + K 4 ) . 



Example Problem 10.8. Consider the system of ODEs: 

f x[(t) =xi-x\ 
x' 2 (t) = t + X\ + x 3 
x' 3 (t) = x 2 — x\ 
*i(0) = 1, 
x 2 (0) = 0, 
I x 3 (0) = 1. 

Approximate X(0.1) by taking a single step of size 0.1, for Euler's Method, and the Runge-Kutta 
Method of order 2. 
Solution: We write 





x 2 (t) - xl(t) 




" 1 " 


F{t,X(t)) = 


t + Xi(t) + x 3 {t) 


X(0) = 







x 2 {t)-x\{t) 




1 



Thus we have 

F(0,X(0)) 

Euler's Method then makes the approximation 

X(0.1) <- X(0) + O.LF(0, X(0)) 
The Runge-Kutta Method computes: 
K x <- 0.1F(0,X(0)) = 



-1 

2 

-1 



" 1 " 




" -0.1 " 




" 0.9 " 





+ 


0.2 




0.2 


1 




-0.1 




0.9 



-0.1 
0.2 
-0.1 



and K 2 ^0.1F(0A,X(0) + K 1 ) 



-0.061 
0.19 
-0.061 
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Then 



X(0.1)^X(0) + -(K 1 + K 2 ) 



' 1 " 




" -0.0805 " 




" 0.9195 " 





+ 


0.195 




0.195 


1 




-0.0805 




0.9195 



10.3.3 It's Only Systems 

The fact we can simply solve systems of ODEs using old methods is good news. It allows us to 
solve higher order ODEs. For example, consider the following problem: given f,a,co,ci, . . . , c n , 
find x(t) such that 



^(t) = /(^(t),x'(f),..,xM)(i)), 

x(a) = cq, x'(a) = ci, x"(a) = c 2 , . . . , x^ n ~ 1 \a) = 



c n -i- 



We can put this in terms of a system by letting 



X q — X ^ 3^ ]_ — «^ 5 X 2 — ^ • • • ^ — 1 — *^ ^ ^ • 



Then the ODE plus these n — 1 equations can be written as 



a4_i(*) = /(*,»()(*), zi(i),x 2 (<),..., x„_i(i)) 

*{,(*)= Xi(t) 

xi(t) = i 2 (<) 

<- 2 (*) = X n -l(t) 

, x (a) = c , xi(a) = ci, x 2 (a) = c 2 , . . . , x n _i(a) = c n _i. 



This is just a system of ODEs that we can solve like any other. 

Note that this trick also works on systems of higher order ODEs. That is, we can transform 
any system of higher order ODEs into a (larger) system of first order ODEs. 

Example Problem 10.9. Rewrite the following system as a system of first order ODEs: 



x'>(t) = (t/y(t))+x>(t)-l 

y (*) = (x'(t)+»(t)) 

I x(0) = 1, x'(0) = 0,2/(0) = -1 



Solution: We let xq = x, x± = x', x 2 = y, then we can rewrite as 



f a^t) = (t/x 2 (t))+x 1 (t) - 1 
4(t) = xi(t) 

X 2(*) = ( Xl {t)+x 2 (t)) 

{ x (Q) = 1, xi(0) = 0,x 2 (0) = 
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10.3.4 It's Only Autonomous Systems 

In fact, we can use a similar trick to get rid of the time-dependence of an ODE. Consider the first 
order system 



fi (t,X 1 ,X 2 , ...,x n ), 
f 2 (t,x 1 ,x 2 ,...,x n ), 
f 3 (t,Xl,X 2 , ...,x n ), 



X n — f n (t,X±, X 2 , . . . , X n ) , 

xi(a) = ci, x 2 (a) = c 2 , . . .,x n (a) = c n . 

We can get rid of the time dependence and thus make the system autonomous by making t 
another variable function to be found. We do this by letting xq = t, then adding the equations 
x'q = 1, and xo(a) = a, to get the system: 



x'l = fl (X ,X!,X 2 , 

x' 2 = f 2 (x ,xi,x 2 , 

x' 3 = h {X ,X!,X 2 , 



X n — f n (xq, X\, X 2 , . . . , Xn) , 

k x (a) = a, xi(a) = ci, x 2 (a) = c 2 , . . .,x n (a) = c n . 
An autonomous system can be written in the more elegant form: 

X' = F(X) 
X (a) = C. 

Moreover, we can consider the phase curve of an autonomous system. 

One can consider X to be the position of a particle which moves through R n over time. For 
an autonomous system of ODEs, the movement of the particle is dependent only on its current 
position and not on time. 1 A phase curve is a flow line in the vector field F. You can think of 
the vector field F as being the velocity, at a given point, of a river, the flow of which is time 
independent. Under this analogy, the phase curve is the path of a vessel placed in the river. In our 
deterministic world, phase curves never cross. This is because at the point of crossing, there would 
have to be two different values of the tangent of the curve, i.e., the vector field would have to have 
two different values at that point. 

The following example illustrates the idea of phase curves for autonomous ODE. 

Example 10.10. Consider the autonomous system of ODEs, without an initial value: 

x[(t) = -2(x 2 -2.4) 
x' 2 (t) = 3(xi - 1.2) 



We can rewrite this ODE as 



X'(t) = F(X(t)) 



1 Note that if you have converted a system of n equations with a time dependence to an autonomous system, then 
the autonomous system should be considered a particle moving through R n+1 . In this case time becomes one of the 
dimensions, and thus we speak of position, instead of time. 
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F(X) = 

This vector field is plotted in Figure 10.6. 



-2 
3 



4.8 
-3.6 




Figure 10.6: The vector field F (X) is shown in R 2 . The tips of the vectors are shown as circles. 
(Due to budget restrictions, arrowheads were not available for this figure.) 



You should verify that this ODE is solved by 

X(t) = 



V2cos(VQt + t ) 





" 1.2 " 


+ 


2.4 



for any r > 0, and to £ (0, 2ir] . Given an initial value, the parameters r and to are uniquely 
determined. Thus the trajectory of X over time is that of an ellipse in the plane. 2 Thus the 
family of such ellipses, taking all r > 0, form the phase curves of this ODE. We would expect the 
approximation of an ODE to never cross a phase curve. This is because the phase curve represents 
the 

Euler's Method and the second order Runge-Kutta Method were used to approximate the ODE 
with initial value 

" 1.5 



X(0) 



1.5 



2 In the case r = 0, the ellipse is the "trivial" ellipse which consists of a point. 
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Figure 10.7: The simulated phase curves of the system of Example 10.10 are shown for (a) Euler's 
Method and (b) the Runge-Kutta Method of order 2. The Runge-Kutta Method used larger step 
size, but was more accurate. The actual solution to the ODE is an ellipse, thus the "spiralling" 
observed for Euler's Method is spurious. 
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Euler's Method was used for step size h = 0.005 for 3000 steps. The Runge-Kutta Method was 
employed with step size h = 0.015 for 3000 steps. The approximations are shown in Figure 10.7. 

Given that the actual solution of this initial value problem is an ellipse, we see that Euler's 
Method performed rather poorly, spiralling out from the ellipse. This must be the case for any step 
size, since Euler's Method steps in the direction tangent to the actual solution; thus every step of 
Euler's Method puts the approximation on an ellipse of larger radius. Smaller stepsize minimizes 
this effect, but at greater computational cost. 

The Runge-Kutta Method performs much better, and for larger step size. At the given resolu- 
tion, no spiralling is observable. Thus the Runge-Kutta Method outperforms Euler's Method, and 
at lesser computational cost. 3 



3 Because the Runge-Kutta Method of order two requires two evaluations of F , whereas Euler's Method requires 
only one, per step, it is expected that they would be 'evenly matched' when Runge-Kutta Method uses a step size 
twice that of Euler's Method. 
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Exercises 

(10.1) Given the ODE: 

r x '(t) = t 2 - x 2 

\ x(l) = 1 

Approximate x(l.l) by using one step of Euler's Method. Approximate x(l.l) using one 
step of the second order Runge-Kutta Method. 

(10.2) Given the ODE: 

x[(t) = X\X 2 + t 
x' 2 (t) = 2X2 — x 2 
xi(l) = 
I x 2 (l) = 3 

Approximate X(l.l) by using one step of Euler's Method. Approximate (1.1) using one 
step of the second order Runge-Kutta Method. 

(10.3) Consider the separable ODE 

x'(t) = x(t)g(t), 

for some function g. Use backward's Euler (or, alternatively, the right endpoint rule for 
approximating integrals) to derive the approximation scheme 

/ , \ x (t) 

J < t + fc) ^ i-MtW 

What could go wrong with using such a scheme? 

(10.4) Consider the ODE 

x '{t) = x{t)+t 2 . 

Using Taylor's Theorem, derive a higher order approximation scheme to find x(t + h) based 
on x(t), t, and h. 

(10.5) Which of the following ODEs can you show are stable? Which are unstable? Which seem 
ambiguous? 

(a) x'(t) = —x — arctanx (b) x'(t) = 2x + tan a; (c) x'(t) = —Ax — e x (d) x'(t) = x + x 3 
(e) x'it) = t + cos x 

(10.6) Rewrite the following higher order ODE as a system of first order ODEs: 

' x"(t) = x(t)- sin £'(*), 

< x(0) = 0, 

k X'(0) = 7T. 

(10.7) Rewrite the system of ODEs as a system of first order ODEs 

' x"(t)=y(t)+t-x(t), 
y'(t) = x'{t) + x(t) - 4, 

< x'(0) = 1, 
x(0) = 2, 

I y(0) = 0. 

(10.8) Rewrite the following higher order system of ODEs as a system of first order ODEs: 

f x"'(t)=y"(t)-x'(t), 

y"'(t)=x"(t)+y'(t), 

X (0) = x"(0) = y'(0) = -1, 
k x'(0) = y(0) = y"(0) = 1. 
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(10.9) Implement Euler's Method for approximating the solution to 

= /(*,*(*)), 
x(a) = c. 

Your m-file should have header line like: 
function xfin = euler (f , a, c ,h, steps) 

where xfin should be the approximate solution to the ODE at time a + steps *h. 

(a) Run your code with f (t, x) = — x, using a = / c, and using h * steps = 1.0. In this 
case the actual solution is 

xfin = ce" 1 . 

This ODE is stable, so your approximation should be good (cf. Example 10.3). Exper- 
iment with different values of h. 

(b) Run your code with f(t,x) = x, using a = / c, and using h * steps = 1.0. In this 
case the actual solution is 

xfin = ce. 

Your actual solution may be a rather poor approximation. Prepare a plot of error 
versus h for varying values of h. (cf. Figure 7.1 and Figure 10.2) 

(c) Run your code with f (t, x) = 1/ (1 — x) , using a = 0, c = 0.1, and using h*steps = 2.0. 
Plot your approximations for varying values of steps. For example, try steps values of 
35, 36, and 90, 91. Also try large values of steps, like 100, 500, 1000. Can you guess what 
the actual solution is supposed to look like? Can you explain the poor performance of 
Euler's Method for this ODE? 

(10.10) Implement the Runge-Kutta Method of order two for approximating the solution to 

= /(*,*(*)), 
x(a) = c. 

Your m-file should have header line like: 
function xfin = rkm(f , a, c ,h, steps) 

where xfin should be the approximate solution to the ODE at time a + steps *h. 
Test your code with the ODEs from the previous problem. Does the Runge-Kutta Method 
perform any better than Euler's Method for the last ODE? 
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Old Exams 



A.l First Midterm, Fall 2003 

PI (7 pnts) Clearly state Taylor's Theorem for f(x + h). Include all hypotheses, and state the 

conditions on the variable that appears in the error term. 
P2 (7 pnts) Write out the Lagrange Polynomial £o{x) for nodes xo,xi, . . . ,x n . That is, write the 

polynomial that has value 1 at xq and has value at x\, X2, ■ ■ ■ , x n . 
P3 (7 pnts) Write the Lagrange form of the polynomial of lowest degree that interpolates the 

following values: 



Xk 


-1 


1 


Vk 


11 


17 



P4 (14 pnts) Use the method of Bisection for finding the zero of the function f(x) = x 3 — 3x 2 + 

5x — |. Let ao = 0, and &o = 2. What are 02, 62? 
P5 (21 pnts) • Write out the iteration step of Newton's Method for finding a zero of the func- 
tion f(x). Your iteration should look like 

x k +i = ?? 

• Write the Newton's Method iteration for finding y/b. 

• Letting xq = 1, find the x\ and X2 for your iteration. Note that y/5 ~ 2.236. 

P6 (7 pnts) Write out the iteration step of the Secant Method for finding a zero of the function 
f(x). Your iteration should look like 

x k+ i =?? 

P7 (21 pnts) Consider the following approximation: 

, ^ /(x + fe)-5/(x) + 4/(x+|) 
3h 

• Derive this approximation using Taylor's Theorem. 

• Assuming that f{x) has bounded derivatives, give the accuracy of the above approxima- 
tion. Your answer should be something like O 

• Let f(x) = x 2 . Approximate /'(0) with this approximation, using h = |. 
P8 (21 pnts) Consider the data: 



k 





1 


2 


Xk 


3 


1 


-1 


f(x k ) 


5 


15 


9 
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• Construct the divided differences table for the data. 

• Construct the Newton form of the polynomial p(x) of lowest degree that interpolates 
f(x) at these nodes. 

• Suppose that these data were generated by the function f(x) = x 3 — 5x 2 + 2x + 17. 
Find a numerical upper bound on the error \p(x) — f{x)\ over the interval [—1,3]. Your 
answer should be a number. 



A. 2 Second Midterm, Fall 2003 



PI (6 pnts) Write the Trapezoid rule approximation for f(x)dx, using the nodes {xi}™ =0 . 
P2 (6 pnts) Suppose Gaussian Elimination with scaled partial pivoting is applied to the matrix 



0.0001 


-0.1 


-0.01 


0.001 


43 


91 


-43 


1 


-12 


26 


-1 





-7 


l 

2 


-10 


7T 



Which row contains the first pivot? You must show your work. 
P3 (14 pnts) Let cf)(n) be some quadrature rule that divides an interval into 2 n equal subintervals. 
Suppose 



f f(x) dx = <j>(n) + a 5 h 5 n + a w h™ + a 15 h£ + ..., 

J a 



where h n = Using evaluations of (/>(•), devise a "Romb erg-style" quadrature rule that is 
a O approximation to the integral. 

P4 (14 pnts) Compute the LU factorization of the matrix 



A = 



1 -1 
2 

2 2-1 



using naive Gaussian Elimination. Show all your work. Your final answer should be two 
matrices, L and U. Verify your answer, i.e., verify that A = LU. 
P5 (14 pnts) Run one iteration of either Jacobi or Gauss Seidel on the system: 



1 


3 


4 " 




' -3 " 


-2 


2 


1 


x = 


-3 


-1 


3 


-2 




1 



Start with = [1 1 1] T . Clearly indicate your answer x^ > . Show your work. You must 
indicate which of the two methods you are using. 
P6 (21 pnts) Derive the two-node Chebyshev Quadrature rule: 



J 1 J{x)dx^c[f(x {i ) + f(x 1 )}. 



Make this rule exact for all polynomials of degree < 2. 
P7 (28 pnts) Consider the quadrature rule: 

f 1 f(x) dx « zof{0) + zi/(l) + z 2 f'(0) + z 3 f'(l). 
Jo 
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• Write four linear equations involving Zi to make this rule exact for polynomials of the 
highest possible degree. 

• Solve your system for z using naive Gaussian Elimination. 



A. 3 Final Exam, Fall 2003 



PI (4 pnts) Clearly state Taylor's Theorem for f(x + h). Include all hypotheses, and state the 

conditions on the variable that appears in the error term. 
P2 (4 pnts) State the conditions that characterize S(x) as a natural cubic spline on the interval 

[a, b]. Let the knots be a = to < t\ < ti < . . . < t n = b. 
P3 (4 pnts) Let xq = 3,x± = — 5,X2 = 0. Write out the Lagrange Polynomial ^o(^) fo r these 

nodes; that is, write the polynomial that has value 1 at xo and has value at xi,X2- 
P4 (4 pnts) Consider the ODE: 

x'(t) = f(t,x(t)) 
x{a) = c 

Write out the Euler Method to step from x(t) to x(t + h). 
P5 (7 pnts) Compute the LU factorization of the matrix 



A = 



-1 1 
1 1 -1 
-2 2 



using naive Gaussian Elimination. Show all your work. Your final answer should be two 
matrices, L and U. Verify your answer, i.e., verify that A = LU. 
P6 (9 pnts) Let a be given. Determine a quadrature rule of the form 



J —t 



f{x) dx « A/(0) + Bf"(-a) + Cf"(a) 



that is exact for polynomials of highest possible degree. What is the highest degree polynomial 
for which this rule is exact? 
P7 (9 pnts) Suppose that 4>{n) is some computational function that approximates a desired quan- 
tity, L. Furthermore suppose that 



4>(n) 



_1 _2 _3 

L + a\n 2 + a2fi 2 + a^n 2 + . . . 



Combine two evaluations of </>(■) in the manner of Richardson to get aO(n x ) approximation 
to L. 

P8 (9 pnts) Find the least-squares best approximate of the form y = In (cx) to the data 



x 


Xl 






2/o 


V\ 




Vn 



Caution: this should not be done with basis vectors. 
P9 (9 pnts) Let 

T = {/(x) = CoX + Cl logX + C2 COSX I Co, C\,C2 6l). 

Set up (but do not attempt to solve) the normal equations to find the least-squares best / G T 
to approximate the data: 



x 


Xl 






2/o 


2/1 




2/n 
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The normal equations should involve the vector c = [co c\ C2] , which uniquely determines a 
function in T . 
P10 (9 pnts) Consider the following approximation: 



/"(*) 



3/(x - fe) - 4/(x) + /(x + 3h) 
6h 2 



• Derive this approximation using Taylor's Theorem. 

• Assuming that /(x) has bounded derivatives, give the accuracy of the above approxima- 
tion. Your answer should be something like O 

Pll (11 pnts) Consider the ODE: 

x'(t) = x{t)g(t) 
x(0) = 1 

Use the Fundamental Theorem of Calculus to write a step update of the form: 



rt+h 

:{t + h) = x(t) + J ?? 



• Approximate the integral using the Trapezoid Rule to get an explicit step update of the 
form 

x(t + h) «- Q(t,x{t),h) 

• Suppose that g(t) > m > for all t. Do you see any problem with this stepping method 
iih> 2/m? 

Hint: the actual solution to the ODE is 



P12 (9 pnts) Consider the ODE: 



x(t) = e& 9 ^ Ar . 



x'{t) = e x ^ +sm(x(t)) 
x(0) = 



Write out the Taylor's series method of order three to approximate x(t + h) given x(t). It is 
permissible to write the update as a small program like: 

At) ^ Qo(x(t)) 

X"(t) <r- Q^xit^At)) 



x(t + h) «- R(x(t),x'{t),x"{t),...) 
rather than write the update as a monstrous equation of the form 

x{t + h) *- M(x(t)). 

P13 (15 pnts) Consider the data: 



k 





1 


2 


Xk 


-0.5 





0.5 


f(Xk) 


5 


15 


9 



• Construct the divided differences table for the data. 
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• Construct the Newton form of the polynomial p(x) of lowest degree that interpolates 
f(x) at these nodes. 

• Suppose that these data were generated by the function f(x) = 40x 3 — 32x 2 — 6x + 15. 
Find a numerical upper bound on the error \p(x) — f(x)\ over the interval [—0.5,0.5]. 
Your answer should be a number. 

P14 (15 pnts) • Write out the iteration step of Newton's Method for finding a zero of the 
function f(x). Your iteration should look like 

Xk+l =?? 

• Write the Newton's Method iteration for finding ^/25/4. 

• Letting xq = |, find the x\ and X2 for your iteration. Note that ^25/4 ~ 1.842. 

P15 (15 pnts) • Let x be some given number, and let /(•) be a given, "black box" function. 

Using Taylor's Theorem, derive a O (h) approximation to the derivative f'{x). Call your 
approximation <fi(h); the approximation should be defined in terms of /(•) evaluated at 
some values, and should probably involve x and h. 

• Let f(x) = x 2 + 1. Approximate /'(0) using your c/)(h), for h=\. 

• For the same f(x) and h, get a O (/i 2 ) approximation to /'(0) using the method of 
Richardson's Extrapolation; your answer should probably include a table like this: 

£>(0,0) 

£(1,0) £(1,1) 

A.4 First Midterm, Fall 2004 

[Definitions] Answer no more than two of the following. 

PI (12 pnts) Clearly state Taylor's Theorem for f(x + h). Include all hypotheses, and state the 

conditions on the variable that appears in the error term. 
P2 (12 pnts) Write the general form of the iterated solution to the problem Ax = b. Let Q be 

your "splitting matrix," and use the factor u. Your solution should look something like: 

??a; (fe) =??a; (fc - 1) +?? or x^ =??a; (fe - 1) +?? 

For one of the following variants, describe what Q is, in relation to A : Richardson's, Jacobi, 
Gauss-Seidel. 

P3 (12 pnts) Write the iteration for Newton's Method or the Secant Method for finding a root to 
the equation f(x) = 0. Your solution should look something like: 

?? 

Xk+l =?? + ^ 

[Problems] Answer no more than four of the following. 

PI (14 pnts) Perform bisection on the function graphed below. Let the first interval be [0,256]. 

Give the second, third, fourth, and fifth intervals. 
P2 (14 pnts) Give a O (/i 4 ) approximation to sin (h + 7r/2) . Find a reasonable C such that the 

truncation error is no greater than Ch 4 . 
P3 (14 pnts) Use Newton's Method to find a good approximation to \/24. Use xq = 3. That is, 

iteratively define a sequence with xq = 3, such that Xk — ► 
P4 (14 pnts) One of the roots of x 2 — bx + c = as found by the quadratic formula is subject 

to subtractive cancellation when \b\ 3> |c| . Which root is it? Rewrite the expression for that 

root to eliminate subtractive cancellation. 




P5 (14 pnts) Find the LU decomposition of the following matrix by Gaussian Elimination with 
naive pivoting. 

"3 2 4 
-6 1 -6 
12 -12 10 _ 

P6 (14 pnts) Consider the linear equation: 



" 6 


1 


1 " 




5 " 


2 


4 





x = 


-6 


_ 1 


2 


6 




3 



Letting x^ = [11 1] T , find a^ 1 ) using either the Gauss-Seidel or Jacobi iterative schemes. 
Use weighting u = 1. You must indicate which scheme you are using. Show all your work. 
[Questions] Answer only one of the following. 

PI (20 pnts) Suppose that the eigenvalues of A are in [a,(3] with < a < (3. Find the to that 
minimizes ]|l — u;A 2 || 2 . What is the minimal value of ||l — <^A 2 || 2 for this optimal uil For full 
credit you need to make some argument why this cj minimizes this norm. 

P2 (20 pnts) Let f(x) = have a unique root r. Recall that the form of the error for Newton's 
Method is 

p _ -/"(&) , 

efc+1 - m^) k ' 

where e& = r — x^. Suppose that f(x) has the properties that f(x) > m > for all x, and 
\f"( x )\ ^ M for all x. Find 5 > such that je^j < 5 insures that x n — ► r. Justify your answer. 
P3 (20 pnts) Suppose that f(r) = = f'(r) / f"(r), and f"(x) is continuous. Letting = x^ — r, 
show that for Newton's Method, 

1 

when |efc| is sufficiently small. (Hint: Use Taylor's Theorem twice.) Is this faster or slower 
convergence than for Newton's Method with a simple root? 



A. 5 Second Midterm, Fall 2004 

[Definitions] Answer no more than two of the following three. 

PI (10 pnts) For a set of distinct nodes xo,xi, . . . ,xi3, what is the Lagrange Polynomial ^(a;)? 
You may write your answer as a product. What degree is this polynomial? 
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P2 (10 pnts) Complete this statement of the polynomial interpolation error theorem: Let p be the 
polynomial of degree at most n interpolating function f at the n+1 distinct nodes Xo,x±, . . . ,x n 
on [a,b]. Let /( ra+1 ) be continuous. Then for each x G [a, b] there is some £ G [a, b] such that 

P3 (10 pnts) State the conditions which define the function S(x) as a natural cubic spline on the 
interval [a, b] . 

[Problems] Answer no more than five of the following six. 

PI (16 pnts) How large must n be to interpolate the function e x to within 0.001 on the interval 

[0, 2] with the polynomial interpolant at n equally spaced nodes? 
P2 (16 pnts) Let 4>(h) be some computable approximation to the quantity L such that 



,15 



4>(h) = L + a 5 h b + ai /i + a i5 h 

Combine evaluations of </>(•) to devise a O (h w ) approximation to L. 
P3 (16 pnts) Derive the constants A and B which make the quadrature rule 



j: 



f(x) dx « Af 



15 



+ Bf(Q) + Af 



T5 



exact for polynomials of degree < 2. Use this quadrature rule to approximate 

2 



f 1 

TT = — 
J-l 1 



+ X 2 



■ dx 



P4 (16 pnts) Find the polynomial of degree < 3 interpolating the data 



X 


-l 





1 


2 


y 


-2 


8 





4 



(Hint: using Lagrange Polynomials will probably take too long.) 
P5 (16 pnts) Find constants a, (3 such that the following is a quadratic spline on [0,4]. 



Q(x) = < 



x 2 - 3x + 4 

a(x - l) 2 + f3(x - 1) + 2 

,2 



^x z + x — 4 

P6 (16 pnts) Consider the following approximation: 

f( x + h)-5f(x) + 4f(x+l) 



< x < 1 

1 < x < 3 
3 < x < 4 



3h 



• Derive this approximation via Taylor's Theorem. 

• Assuming that f(x) has bounded derivatives, give the accuracy of the above approxima- 
tion. Your answer should be something like O (h 7 ). 

• Let f(x) = x 2 . Approximate f'(0) with this approximation, using h = |. 
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A.6 Final Exam, Fall 2004 



[Definitions] Answer all of the following. 

PI (10 pnts) Clearly state Taylor's Theorem for f(x + h). Include all hypotheses, and state the 

conditions on the variable that appears in the error term. 
P2 (10 pnts) Write the iteration for Newton's Method or the Secant Method for finding a root to 

the equation f(x) = 0. Your solution should look something like: 



Xk+l =?? + 



?? 
77 



P3 (10 pnts) Write the general form of the iterated solution to the problem Ax = b. Let Q be 
your "splitting matrix," and use the factor ui. Your solution should look something like: 



77 x ( k ) =??a;( A;_1 )+?? 



or 



x (k) = ?? iC ( fc - 1 )_)_?? 



For one of the following variants, describe what Q is, in relation to A : Richardson's, Jacobi, 
Gauss-Seidel. 

P4 (10 pnts) Define what it means for the function f(x) € T to be the least squares best approx- 
imant to the data 



Xq 


Xi 






yo 


yi 




Vn 



[Problems] Answer all of the following. 

PI (10 pnts) Rewrite this system of higher order ODEs as a system of first order ODEs: 



x"(t) 


= x'{t) [x(t) + ty{t)\ 


y"(t) 


= y'(t) \y{t) - tx(t)} 


s'(0) 


= 1 


< x(0) = 


= 4 


1/(0) 


= -3 


y(o) -- 


= -2 



P2 (10 pnts) Consider the ODE: 



V(t) =t(x + l) 2 
x(2) = 1/2 



Use Euler's Method to find an approximate value of x(2.1) using a single step. 
P3 (10 pnts) Let (j>{h) be some computable approximation to the quantity L such that 



4>{h) = L + a 4 /i 4 + a 8 /i 8 + a 12 h lz + ... 

Combine evaluations of <p(-) to devise a O (/i 8 ) approximation to L. 
P4 (15 pnts) Consider the approximation 



/'(*) 



9 f{x + h)- 8f(x) - f{x - 3h) 
12h 



(a) Derive this approximation using Taylor's Theorem. 

(b) Assuming that f{x) has bounded derivatives, give the accuracy of the above approxima- 
tion. Your answer should be something like O (h 7 ). 
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P5 (15 pnts) Find constants, a, (3, 7 such that the following is a natural cubic spline on [1,3]. 

Q(x) 



'a (x - l) 3 + (x - l) 2 + 12 (x - 1) + 1 : 1 < x < 2 



7 x d + 18x^ - 30x + 19 



: 2 < x < 3 



P6 (15 pnts) Devise a subroutine that, given an input number z, calculates an approximate to 

via Newton's Method. Write your subroutine in pseudo code, Matlab, or C/C++. Include 
reasonable termination conditions so your subroutine will terminate if it finds a good approx- 
imate solution. Your subroutine should use only the basic arithmetic operations of addition, 
subtraction, multiplication, division; in particular your subroutine should not have access to 
a cubed root or logarithm operator. 

P7 (15 pnts) Consider the ODE: 

jx'(t) = cos (x(t)) + sin (x(t)) 
\x(0) = 1 

Use Taylor's theorem to devise a O (/i 3 ) approximant to x(t + h) which is a function of x, t, 
and h only. You may write your answer as a program: 

x >(t) <- Q (x(t)) 
a/'(t)<-Q 1 (x(t),x'(t)) 

x <- R(x{t),x'{t),x"(t),...) 

such that x(t + h) = x + O (h 3 ) . 
[Questions] Answer all of the following. 
PI (20 pnts) Consider the "quadrature" rule: 

/' f(x) dx « A f(0) + Ai/(1) + A 2 f'(l) + A 3 f"(l) 
Jo 

(a) Write four linear equations involving A^ to make this rule exact for polynomials of the 
highest possible degree. 

(b) Find the constants Ai using Gaussian Elimination with naive pivoting. 
P2 (20 pnts) Consider the data: 



k 





1 


2 


x k 





0.25 


0.5 


f(x k ) 


1.5 


1 


0.5 



(a) Construct the divided differences table for the data. 

(b) Construct the Newton form of the polynomial p(x) of lowest degree that interpolates 
/(x) at these nodes. 

(c) What degree is your polynomial? Is this strange? 

(d) Suppose that these data were generated by the function 



f(x) = 1 + 



COS 2"7TX 



Find an upper bound on the error \p(x) — f(x)\ over the interval [0,0.5]. Your answer 
should be a number. 
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P3 (20 pnts) Let 

T = {f(x) = coVx + ci£(x) + c 2 e x | c , c±, c 2 eK}, 

where £(x) is some black box function. We are concerned with finding the least-squares best 
/ G T to approximate the data 



x 


Xl 




Xn 


yo 


yi 




Vn 



(a) Set up (but do not attempt to solve) the normal equations to find the the vector c = 
[co ci C2] T , which determines the least squares approximant in T . 

(b) Now suppose that the function l(x) happens to be the Lagrange Polynomial associated 
with X7, among the values XQ,X\,...,x n . That is, £(xj) = 5jj. Prove that f(xj) = j/7, 
where / is the least squares best approximant to the data. 

P4 (10 pnts) State some substantive question which you thought might appear on this exam, but 
did not. Answer this question (correctly). 
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GNU Free Documentation License 



Version 1.2, November 2002 
Copyright ©2000,2001,2002 Free Software Foundation, Inc. 

59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 

Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed. 

Preamble 

The purpose of this License is to make a manual, textbook, or other functional and useful document "free" in the sense of 
freedom: to assure everyone the effective freedom to copy and redistribute it, with or without modifying it, either commercially 
or noncommercially. Secondarily, this License preserves for the author and publisher a way to get credit for their work, while 
not being considered responsible for modifications made by others. 

This License is a kind of " copy left" , which means that derivative works of the document must themselves be free in the 
same sense. It complements the GNU General Public License, which is a copyleft license designed for free software. 

We have designed this License in order to use it for manuals for free software, because free software needs free documentation: 
a free program should come with manuals providing the same freedoms that the software does. But this License is not limited 
to software manuals; it can be used for any textual work, regardless of subject matter or whether it is published as a printed 
book. We recommend this License principally for works whose purpose is instruction or reference. 

1. APPLICABILITY AND DEFINITIONS 

This License applies to any manual or other work, in any medium, that contains a notice placed by the copyright holder 
saying it can be distributed under the terms of this License. Such a notice grants a world-wide, royalty-free license, unlimited 
in duration, to use that work under the conditions stated herein. The "Document", below, refers to any such manual or work. 
Any member of the public is a licensee, and is addressed as "you". You accept the license if you copy, modify or distribute the 
work in a way requiring permission under copyright law. 

A "Modified Version" of the Document means any work containing the Document or a portion of it, cither copied 
verbatim, or with modifications and/or translated into another language. 

A "Secondary Section" is a named appendix or a front-matter section of the Document that deals exclusively with the 
relationship of the publishers or authors of the Document to the Document's overall subject (or to related matters) and contains 
nothing that could fall directly within that overall subject. (Thus, if the Document is in part a textbook of mathematics, a 
Secondary Section may not explain any mathematics.) The relationship could be a matter of historical connection with the 
subject or with related matters, or of legal, commercial, philosophical, ethical or political position regarding them. 

The "Invariant Sections" are certain Secondary Sections whose titles are designated, as being those of Invariant Sections, 
in the notice that says that the Document is released under this License. If a section does not fit the above definition of Secondary 
then it is not allowed to be designated as Invariant. The Document may contain zero Invariant Sections. If the Document does 
not identify any Invariant Sections then there are none. 

The "Cover Texts" are certain short passages of text that are listed, as Front-Cover Texts or Back-Cover Texts, in 
the notice that says that the Document is released under this License. A Front-Cover Text may be at most 5 words, and a 
Back-Cover Text may be at most 25 words. 

A "Transparent" copy of the Document means a machine-readable copy, represented in a format whose specification is 
available to the general public, that is suitable for revising the document straightforwardly with generic text editors or (for 
images composed of pixels) generic paint programs or (for drawings) some widely available drawing editor, and that is suitable 
for input to text formatters or for automatic translation to a variety of formats suitable for input to text formatters. A copy 
made in an otherwise Transparent file format whose markup, or absence of markup, has been arranged to thwart or discourage 
subsequent modification by readers is not Transparent. An image format is not Transparent if used for any substantial amount 
of text. A copy that is not "Transparent" is called "Opaque". 
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Examples of suitable formats for Transparent copies include plain ASCII without markup, Texinfo input format, LaTeX 
input format, SGML or XML using a publicly available DTD, and standard-conforming simple HTML, PostScript or PDF 
designed for human modification. Examples of transparent image formats include PNG, XCF and JPG. Opaque formats 
include proprietary formats that can be read and edited only by proprietary word processors, SGML or XML for which the 
DTD and/or processing tools are not generally available, and the machine-generated HTML, PostScript or PDF produced by 
some word processors for output purposes only. 

The "Title Page" means, for a printed book, the title page itself, plus such following pages as are needed to hold, legibly, 
the material this License requires to appear in the title page. For works in formats which do not have any title page as such, 
"Title Page" means the text near the most prominent appearance of the work's title, preceding the beginning of the body of 
the text. 

A section "Entitled XYZ" means a named subunit of the Document whose title either is precisely XYZ or contains XYZ 
in parentheses following text that translates XYZ in another language. (Here XYZ stands for a specific section name mentioned 
below, such as "Acknowledgements", "Dedications", "Endorsements", or "History".) To "Preserve the Title" of 
such a section when you modify the Document means that it remains a section "Entitled XYZ" according to this definition. 

The Document may include Warranty Disclaimers next to the notice which states that this License applies to the Document. 
These Warranty Disclaimers are considered to be included by reference in this License, but only as regards disclaiming warranties: 
any other implication that these Warranty Disclaimers may have is void and has no effect on the meaning of this License. 

2. VERBATIM COPYING 

You may copy and distribute the Document in any medium, cither commercially or noncommercially, provided that this 
License, the copyright notices, and the license notice saying this License applies to the Document are reproduced in all copies, 
and that you add no other conditions whatsoever to those of this License. You may not use technical measures to obstruct or 
control the reading or further copying of the copies you make or distribute. However, you may accept compensation in exchange 
for copies. If you distribute a large enough number of copies you must also follow the conditions in section 3. 

You may also lend copies, under the same conditions stated above, and you may publicly display copies. 

3. COPYING IN QUANTITY 

If you publish printed copies (or copies in media that commonly have printed covers) of the Document, numbering more 
than 100, and the Document's license notice requires Cover Texts, you must enclose the copies in covers that carry, clearly and 
legibly, all these Cover Texts: Front-Cover Texts on the front cover, and Back-Cover Texts on the back cover. Both covers must 
also clearly and legibly identify you as the publisher of these copies. The front cover must present the full title with all words of 
the title equally prominent and visible. You may add other material on the covers in addition. Copying with changes limited to 
the covers, as long as they preserve the title of the Document and satisfy these conditions, can be treated as verbatim copying 
in other respects. 

If the required texts for either cover are too voluminous to fit legibly, you should put the first ones listed (as many as fit 
reasonably) on the actual cover, and continue the rest onto adjacent pages. 

If you publish or distribute Opaque copies of the Document numbering more than 100, you must either include a machine- 
readable Transparent copy along with each Opaque copy, or state in or with each Opaque copy a computer-network location from 
which the general network-using public has access to download using public-standard network protocols a complete Transparent 
copy of the Document, free of added material. If you use the latter option, you must take reasonably prudent steps, when you 
begin distribution of Opaque copies in quantity, to ensure that this Transparent copy will remain thus accessible at the stated 
location until at least one year after the last time you distribute an Opaque copy (directly or through your agents or retailers) 
of that edition to the public. 

It is requested, but not required, that you contact the authors of the Document well before redistributing any large number 
of copies, to give them a chance to provide you with an updated version of the Document. 

4. MODIFICATIONS 

You may copy and distribute a Modified Version of the Document under the conditions of sections 2 and 3 above, provided 
that you release the Modified Version under precisely this License, with the Modified Version filling the role of the Document, 
thus licensing distribution and modification of the Modified Version to whoever possesses a copy of it. In addition, you must 
do these things in the Modified Version: 

A. Use in the Title Page (and on the covers, if any) a title distinct from that of the Document, and from those of previous 
versions (which should, if there were any, be listed in the History section of the Document). You may use the same title 
as a previous version if the original publisher of that version gives permission. 

B. List on the Title Page, as authors, one or more persons or entities responsible for authorship of the modifications in the 
Modified Version, together with at least five of the principal authors of the Document (all of its principal authors, if it 
has fewer than five), unless they release you from this requirement. 

C. State on the Title page the name of the publisher of the Modified Version, as the publisher. 

D. Preserve all the copyright notices of the Document. 

E. Add an appropriate copyright notice for your modifications adjacent to the other copyright notices. 
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F. Include, immediately after the copyright notices, a license notice giving the public permission to use the Modified Version 
under the terms of this License, in the form shown in the Addendum below. 

G. Preserve in that license notice the full lists of Invariant Sections and required Cover Texts given in the Document's 
license notice. 

H. Include an unaltered copy of this License. 

I. Preserve the section Entitled "History", Preserve its Title, and add to it an item stating at least the title, year, new 
authors, and publisher of the Modified Version as given on the Title Page. If there is no section Entitled "History" in 
the Document, create one stating the title, year, authors, and publisher of the Document as given on its Title Page, 
then add an item describing the Modified Version as stated in the previous sentence. 

J. Preserve the network location, if any, given in the Document for public access to a Transparent copy of the Document, 
and likewise the network locations given in the Document for previous versions it was based on. These may be placed 
in the "History" section. You may omit a network location for a work that was published at least four years before the 
Document itself, or if the original publisher of the version it refers to gives permission. 

K. For any section Entitled "Acknowledgements" or "Dedications", Preserve the Title of the section, and preserve in the 
section all the substance and tone of each of the contributor acknowledgements and/or dedications given therein. 

L. Preserve all the Invariant Sections of the Document, unaltered in their text and in their titles. Section numbers or the 
equivalent are not considered part of the section titles. 

M. Delete any section Entitled "Endorsements". Such a section may not be included in the Modified Version. 

N. Do not retitle any existing section to be Entitled "Endorsements" or to conflict in title with any Invariant Section. 

O. Preserve any Warranty Disclaimers. 

If the Modified Version includes new front-matter sections or appendices that qualify as Secondary Sections and contain 
no material copied from the Document, you may at your option designate some or all of these sections as invariant. To do this, 
add their titles to the list of Invariant Sections in the Modified Version's license notice. These titles must be distinct from any 
other section titles. 

You may add a section Entitled " Endorsements" , provided it contains nothing but endorsements of your Modified Version by 
various parties-for example, statements of peer review or that the text has been approved by an organization as the authoritative 
definition of a standard. 

You may add a passage of up to five words as a Front-Cover Text, and a passage of up to 25 words as a Back-Cover Text, 
to the end of the list of Cover Texts in the Modified Version. Only one passage of Front-Cover Text and one of Back-Cover 
Text may be added by (or through arrangements made by) any one entity. If the Document already includes a cover text for 
the same cover, previously added by you or by arrangement made by the same entity you are acting on behalf of, you may not 
add another; but you may replace the old one, on explicit permission from the previous publisher that added the old one. 

The author(s) and publisher(s) of the Document do not by this License give permission to use their names for publicity 
for or to assert or imply endorsement of any Modified Version. 

5. COMBINING DOCUMENTS 

You may combine the Document with other documents released under this License, under the terms defined in section 4 
above for modified versions, provided that you include in the combination all of the Invariant Sections of all of the original 
documents, unmodified, and list them all as Invariant Sections of your combined work in its license notice, and that you preserve 
all their Warranty Disclaimers. 

The combined work need only contain one copy of this License, and multiple identical Invariant Sections may be replaced 
with a single copy. If there arc multiple Invariant Sections with the same name but different contents, make the title of each 
such section unique by adding at the end of it, in parentheses, the name of the original author or publisher of that section if 
known, or else a unique number. Make the same adjustment to the section titles in the list of Invariant Sections in the license 
notice of the combined work. 

In the combination, you must combine any sections Entitled "History" in the various original documents, forming one section 
Entitled " History" ; likewise combine any sections Entitled " Acknowledgements" , and any sections Entitled " Dedications" . You 
must delete all sections Entitled "Endorsements". 

6. COLLECTIONS OF DOCUMENTS 

You may make a collection consisting of the Document and other documents released under this License, and replace the 
individual copies of this License in the various documents with a single copy that is included in the collection, provided that 
you follow the rules of this License for verbatim copying of each of the documents in all other respects. 

You may extract a single document from such a collection, and distribute it individually under this License, provided you 
insert a copy of this License into the extracted document, and follow this License in all other respects regarding verbatim 
copying of that document. 



7. AGGREGATION WITH INDEPENDENT WORKS 
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A compilation of the Document or its derivatives with other separate and independent documents or works, in or on a 
volume of a storage or distribution medium, is called an "aggregate" if the copyright resulting from the compilation is not used 
to limit the legal rights of the compilation's users beyond what the individual works permit. When the Document is included 
in an aggregate, this License does not apply to the other works in the aggregate which are not themselves derivative works of 
the Document. 

If the Cover Text requirement of section 3 is applicable to these copies of the Document, then if the Document is less than 
one half of the entire aggregate, the Document's Cover Texts may be placed on covers that bracket the Document within the 
aggregate, or the electronic equivalent of covers if the Document is in electronic form. Otherwise they must appear on printed 
covers that bracket the whole aggregate. 

8. TRANSLATION 

Translation is considered a kind of modification, so you may distribute translations of the Document under the terms of 
section 4. Replacing Invariant Sections with translations requires special permission from their copyright holders, but you may 
include translations of some or all Invariant Sections in addition to the original versions of these Invariant Sections. You may 
include a translation of this License, and all the license notices in the Document, and any Warranty Disclaimers, provided that 
you also include the original English version of this License and the original versions of those notices and disclaimers. In case 
of a disagreement between the translation and the original version of this License or a notice or disclaimer, the original version 
will prevail. 

If a section in the Document is Entitled "Acknowledgements", "Dedications", or "History", the requirement (section 4) to 
Preserve its Title (section 1) will typically require changing the actual title. 

9. TERMINATION 

You may not copy, modify, sublicense, or distribute the Document except as expressly provided for under this License. 
Any other attempt to copy, modify, sublicense or distribute the Document is void, and will automatically terminate your rights 
under this License. However, parties who have received copies, or rights, from you under this License will not have their licenses 
terminated so long as such parties remain in full compliance. 

10. FUTURE REVISIONS OF THIS LICENSE 

The Free Software Foundation may publish new, revised versions of the GNU Free Documentation License from time to 
time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or 
concerns. See http://www.gnu.org/copylcft/. 

Each version of the License is given a distinguishing version number. If the Document specifies that a particular numbered 
version of this License "or any later version" applies to it, you have the option of following the terms and conditions either of 
that specified version or of any later version that has been published (not as a draft) by the Free Software Foundation. If the 
Document does not specify a version number of this License, you may choose any version ever published (not as a draft) by the 
Free Software Foundation. 

ADDENDUM: How to use this License for your documents 

To use this License in a document you have written, include a copy of the License in the document and put the following 
copyright and license notices just after the title page: 

Copyright ©YEAR YOUR NAME. Permission is granted to copy, distribute and/or modify this document 
under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free 
Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of 
the license is included in the section entitled " GNU Free Documentation License" . 



If you have Invariant Sections, Front-Cover Texts and Back-Cover Texts, replace the "with. ..Texts." line with this: 

with the Invariant Sections being LIST THEIR TITLES, with the Front-Cover Texts being LIST, and with the 
Back-Cover Texts being LIST. 

If you have Invariant Sections without Cover Texts, or some other combination of the three, merge those two alternatives 
to suit the situation. 

If your document contains nontrivial examples of program code, we recommend releasing these examples in parallel under 
your choice of free software license, such as the GNU General Public License, to permit their use in free software. 
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Riemann Integrable, 98 
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